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PREFACE 


The NASA Aircraft Noise Prediction Program (ANDPP) was developed 
originally for the prediction of airport community noise from turbofan- 
powered aircraft. The theoretical manuals for the original prediction 
system were published as NASA TM-83199, parts 1 and 2, in February 1982. 
Part 1 describes program modules which define the atmosphere, the flight 
of the aircraft, the propagation of the noise, and the subjective 
effects of the noise on the receiver. Part 2 describes program modules 
which define the turbofan engine noise sources and the airframe noise 
sources of CTOL aircraft. 

The purpose of part 3 of the theoretical manual is to describe 
those additional program modules which are used to define propeller 
noise sources and the propagation of pure tones. Even though this part 
begins with chapter 10 to follow the numbered chapters of parts 1 and 2, 
the manual is written such that the chapters have minimal interdepen- 
dence. The program user may rely on this part of the manual to define 
the propeller noise sources and then refer to parts 1 and 2 to find the 
effect of these sources on the airport community. 
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10. PROPELLER AERODYNAMICS 


10.1 PROPELLER ANALYSIS SYSTEM 


William E. Zorumski 
Langley Research Center 


SUMMARY 

The NASA Aircraft Noise Prediction Program (ANOPP) Propeller 
Analysis System is a set of computational modules for predicting the 
aerodynamics, performance, and noise of propellers. Propeller blade 
geometry is given in terms of blade surface coordinates derived from a 
Joukowski transform of the blade sections. Potential flow around the 
blade sections is computed by Theodorsen's method by using the Kutta 
condition to fix the circulation. Blade boundary layers are computed 
by using the Holstei n-Bohlen method in the laminar region and the 
Truckenbrodt method in the turbulent region. Profile drag is predicted 
by the method of Young and Squires. Performance and induced flow are 
computed by Lock's method with the Prandtl ci rcul ati on ■ function near 
the blade tip. Discrete tone noise is predicted from blade shape and 
aerodynamic loads by using Farassat's methods: the blade surface 

integral method for subsonic propellers and the collapsing sphere 
method for transonic propellers. Broadband trai 1 i ng-edge noise is 
computed by the method of Schlinker and Amiet. The results of this 
prediction system are compared with measurements on two propellers: 
one subsonic and one transonic. Near-field levels on the subsonic 
propeller are accurately predicted if the predicted power coefficient 
is adjusted to match the measured power coefficient of the propeller. 
The lower frequency harmonics of the subsonic propeller spectrum match 
the measured values but the high frequency harmonics are under- 
predicted. This underpredi ction is believed to be caused by the 
omission of unsteady loading effects in the predictions. The far-field 
or flyover noise of the subsonic propeller is scattered by atmospheric 
and ground effects but the general trend of the data indicates over- 
prediction of far-field levels. Transonic propeller noise measured on 
the fuselage of the aircraft is significantly influenced by the refrac- 
tion effects of the boundary layer on the fuselage. When these effects 
are included and when the power is matched, the transonic predictions 
agree with the data except in a small region just behind the propeller 
plane on the aircraft surface. It is believed that fuselage scattering 
effects may be the cause of this discrepancy. 


INTRODUCTION 

Propeller noise prediction is based on two disciplines: aero- 

dynamics and acoustics. NASA has developed the Aircraft Noise 
Prediction Program (ANOPP) for aircraft noise prediction. This part of 
the manual develops the methods for propeller noise prediction by using 
the assumption that these disciplines are separable; that is, the flow 
field can be separated into an aerodynamic part and acoustic part. 

This separation allows the computations to be made sequentially. 
Classical aerodynamic theory is used to find the surface pressures and 
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frictional stresses on the blade surfaces, and then acoustic theories 
are used to predict the noise. 

The ANOPP propeller noise prediction system is diagrammed in 
figure 1. The different computational tasks are assigned to 
independent blocks of computer code called functional modules. These 
modules are managed by an executive system, the ANOPP executive program 
(ref. 1). One group of modules deals with the aerodynamic computa- 
tions. Not shown are modules which compute atmospheric properties and 
flight dynamics. These modules are described in reference 2. The end 
products of the aerodynamics modules are the propeller blade motions 
and loads. Motions include only the aircraft motion and rotational 
effects. Flexing and vibration, while possibly important, are not 
included at this time. The loads are the pressure and frictional 
stress on the propeller blade surface. These loads are generally a 
function of both surface position and time. Given the blade motions 
and loads, it is theoretically possible to predict the noise. As a 
practical matter, however, this prediction is impossible. The loads 
are really nonstationary random processes, and there is presently no 
feasible computational procedure which will produce a complete descrip- 
tion of the noise. The approximation is made that the noise may be 
divided into two parts: discrete and random. The discrete (tone) 
noise is computed directly from the blade motions and loads, whereas 
the randan (broadband) noise is estimated by semi -empi ri cal methods. 

The discrete noise computations are made entirely with time domain 
techniques, and randan noise estimates use a blend of frequency and 
time domain concepts. 

Noise predictions are needed in both the near and far fields. 
Near-field noise predictions are needed to find the noise transmitted 
to the interior of an aircraft. There are two important near-field 
effects on propel 1 er-generated noise. These effects are the refraction 
caused by nonuniform flow over the aircraft and the scattering by the 
aircraft body. Far-field effects are needed to accurately predict 
community noise. These effects are the atmospheric attenuation and 
refraction, and the ground effects of reflection and absorption. 

Further information on the far-field propagation modules may be found 
in reference 2. Modules for the near-field effects are described in 
detail herein. 

Agreement with experiment is the ultimate goal of a prediction 
system. Any prediction method - a guess, a curve fit, or a solution to 
a partial differential equation - is acceptable if it agrees with 
experimental data according to some objective rule. The question is 
how many experiments and how good must the agreement be to prove that a 
prediction is correct. This question raises another question: prove 

it to whom? You cannot prove to a member of the Flat Earth Society 
that the Earth is round. 

Despite these minor difficulties, we gather data from experiments 
in the hope of proving that the prediction system is correct. These 
are described as flight data and tunnel data to denote the two types of 
facilities most frequently used to conduct the experiments. The flight 
data are usually but not always full-scale. Results will be shown 
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later for two propellers: one subsonic design and one transonic 

design. The subsonic propeller is made by the TRW Hartzell Propeller 
Products Division and its noise was measured by a wing-mounted boom 
microphone in flight. The transonic propeller is the NASA SR -3 Propfan 
design. This roughly quarter-seal e propeller was tested in flight 
mounted atop a JetStar aircraft. 
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SYMBOLS 

radius of perfect circle 
Joukowski transform parameter 
blade section chord 
energy dissipation coefficient 

blade section drag coefficient, re section dynamic 
pressure and chord 

skin friction coefficient, re section dynamic pressure 

blade section lift coefficient, re section dynamic 
pressure and chord 

Theodorsen transform coefficients 

power coefficient 

coefficient of pressure, re section dynamic pressure 

wall shear stress coefficient, re local dynamic 
pressure 

ambient speed of sound 
drag 

Prandtl tip vorticity function 
al ti tude 

boundary-1 ayer shape factors, i * j = 1, 2, 3 
Hankel functions of the first and second kinds 
advance ratio, V^/nD 
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wave numbers 
1 i ft 
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loading intensity vector 

component of 1 in radiation direction 

empirical broadband correlation length 

blade section Mach number 

forward Mach number 

tip helical Mach number 

normal Mach number 

radiation Mach number 

tip rotational Mach number 

free-stream or forward Mach number 

ci rcumferenti al harmonic number 

number of blades 

revolutions per second 

normal to blade surface 

overall sound pressure level 

acoustic pressure 

root-mean-square pressure 

torque 

propeller disk radius 
cylindrical coordinates 
radiation vector, x - y 
unit radiation vector 
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sound pressure level at a single frequency 
surface area 

empirical broadband spectrum 

thrust, period of signal, or transfer function 

time 

tangents to blade surface in spanwise and chordwise 
directions, respectively 

velocity at edge of boundary layer 

local Mach number in boundary layer 

free-stream velocity 

normal acoustic velocity 

acoustic velocity vector 

normal acoustic velocity at edge of boundary layer 

complex velocity field around perfect circle 

distance along chord or distance along airfoil surface 
measured from stagnation point; coordinate in flow 
di recti on 

observer position in media-fixed reference frame 

coordinate normal to flow direction 

source position in media-fixed reference frame 

complex plane of the airfoil section 

blade section angle of attack 

angle-of-attack perturbation 

blade pitch at 3/4 span 

circulation around airfoil section or curve of 

intersection of collapsing sphere and propeller blade 
surface 

boundary-layer thickness 

boundary-layer displacement thickness 

boundary-layer momentum thickness 

boundary-layer energy thickness 
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Subscri pts : 
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n 2 > n 3 


?2 


polar coordinate angle difference between near circle 
and perfect circle 

complex plane of airfoil section represented as a 
perfect circle 

complex plane of airfoil section represented as a 
near-ci rcl e 

pitch axis coordinates fixed to propeller; n 2 is 
along pitch axis and n 3 is along forward shaft axis 

angle between normal vector and radiation vector 

advance ratio, M /M. 

co t 

blade surface coordinates; ^ varies with span, c 2 
varies with chord 

ambient density 

propeller solidity 

retarded time 

blade section inflow angle 
elliptic blade thickness function 
propeller angular velocity 
angular frequency 

i nduced 
leading edge 
refraction effect 
stagnation point 
scattering effect 
trailing edge 
transition point 
free stream 
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Notation : 


time dependent 
_ mean 

->■ vector 

1 derivative with respect to argument 


AERODYNAMIC PREDICTION 


Geometry 


Although geometry is given rather than predicted, it is worthwhile 
to carefully consider the way geometric data are utilized. Two 
coordinates are required to specify a point on a surface. It is 
important to choose surface coordinates such that the functions of 
these coordinates will be single valued and free of singularities at 
least up to their second derivatives; that is, functions of the surface 
coordinates will be of class C 2 . The selected surface coordinates 
should also provide a convenient computational grid. 


Figure 2 shows the basic coordinate system (yi,y 2 ,y 3 ) used to 
describe the propeller. At time t = 0, a rotating coordinate system 
( rii , 02 >^ 3 ) coincides with the (y!,y 2 ,y 3 ) system. The shaft axis is 
133 , the blade pitch axis is n 2 , and the m-axis completes an 
orthogonal triad. Section A-A through the blade at constant span 
positions gives the first surface coordinate 


£1 " ^2 

At each section, a Joukowski transform (see ref. 3, for example) 


( 1 ) 


z = ?' 



( 2 ) 


where 


z = m + in 3 (3) 


c 1 = be^ 2 * 1 ^ (4) 

is used to introduce the second surface coordinate ? 2 . The blade 
section elliptical coordinate ip(c 2 ) resulting from this transforma- 
tion is shown in figure 3. The blade section is described by the func- 
tion t|>( 5 2 ). The blade surface is unwrapped by the Joukowski transform 
if the first surface coordinate is chosen to equal n 2 . The 

surface function i|/(?i,? 2 ) is single valued, continuous, and slowly 
varying in the surface coordinates ( 5 1 ,c 2 ) as shown in figure 3. The 
computation grid stretches the region near the leading edge, £ 2 = 0.5, 
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so that aerodynamic functions such as the coefficient of pressure Cp 
will be slowly varying in this region. In this computational space, 
bicubic splines are suitable interpolating functions and will be used 
in all subsequent computations. 


Potential Field 


The potential flow field around each airfoil section 
= Constant is given by a conformal transformation of the flow 
around a perfect cylinder (ref. 3). The blade geometry analysis has 
already produced part of the desired transformati on by mapping the 
airfoil section in the z-plane as shown in figure 4. Recall that the 
airfoil thickness function ^(c 1 ,? 2 ) was generated by inverting the 
Joukowski transformation: 



This transformation may now be used directly to map a given flow around 
the near-circle into a flow around the airfoil section. 


Theodorsen's transformation (ref. 4) maps the ?-plane of the 
perfect circle into the c'-plane of the near-circle 

00 

c ' = 5exp (£)) (6) 

The constants C n in Theodorsen's transformation are found from the 
shape of the airfoil in the c'-plane. After numerically solving for 
these constants, it is found that the trai 1 i ng-edge point of the 
airfoil is displaced by a small angle e T from the real axis of 
the c-plane. This point is required to be a stagnation point for the 
flow around the cylinder in order to satisfy the Kutta condition that 
trai 1 i ng-edge velocities are finite. 


The complex flow function around the perfect circle is 


U(c) 



it- 

^2 / ^ 


(7) 


where M is the local Mach number (for the section), a is the angle 
of attack, a is the radius of the circle, and r is the circula- 
tion. The trai 1 i ng-edge stagnation point is 


1e T 

? ST ae 

which, with W(? ST ) = 0, gives the circulation r as 

r = 4iraM si n( a - e-p) (8) 
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The leading-edge stagnation point, a second solution to the equation 
for WU $ ) is 


C SL 


i ( 2oi-ej+ir) 

ae 


( 9 ) 


The cross product of the velocity vector and the circulation vector 
gi ves the 1 i ft 

C^= 2tt sin (a - e T ) ( 10) 


The coefficient of pressure is found from the velocity in the 
z-plane of the airfoil: 


C 


P 


= i _ 


I w ( z ) | 2 

I T - I 


(ID 


The complex velocity in the z-plane is found by using the derivatives 
of the transformations 


W(z) = 


d? W(e) 

d’? T dz 

dT 


(c' * (12) 


At the stagnation points in the ^-plane, the limit of equation (12) is 
used: 


W(z) = 


d? 

dF 


Lim 

C'< 


W U) 

dz 

dT 


(13) 


Since the Joukowski transformation has a singularity at the trailing 
edge, W(^) will have a finite nonzero limit at this point. This 
limit causes the coefficient of pressure at the trailing edge to have 
magnitude less than 1. 

Figure 5 shows the blade section data computed by the blade 
section potential flow analysis. The coefficient of pressure C p 
has a maximum of unity at the leading-edge stagnation point 5 2 $. 

The lower surface C p is generally positive, decreasing to a small 
value near the trailing edge ? 2 = 1.0, where the airfoil surface 
velocity is near the free-stream velocity. The upper surface C p 
may become negative with high velocities around the highly curved 
leading edge. The upper surface C p then approaches a small 
negative value at the trailing edge £ 2 = 0. The lift coefficient 
Cj_ has a slope of approximately 2-rr when plotted as a function of 
a. 


The function of the Blade Section Aerodynamics Module is 
summarized as follows. Input is the blade shape function ip(z 1 ,Z 2 ) 
and the parameters a and M. The module maps the flow around a 
perfect cylinder into the flow around the blade by using the Theodorsen 
and Joukowski transforms. Outputs are the section lift coefficient C 
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and leading-edge stagnation point c 2 S as a function of span 
position ?! and the parameters a and M. The output coefficient of 
pressure depends on surface position (5 ls 5 2 ) arKl the parameters a 
and M. The stagnation-point location and the coefficient of pressure 
are used in the boundary-1 ayer analysis which follows. 


Boundary Layer 

The Blade Section Boundary-Layer Module computes the 
two-dimensional boundary layer on each airfoil section as shown in 
figure 6. The blade arc length x is measured from the stagnation 
point on the leading edge. The initial portion of the boundary layer 
is laminar. Transition to turbulent flow occurs near the point where 
the external velocity U(x) is a maximum. This turbulent layer 
continues to the trailing edge unless separation occurs. 


The governing equations for the boundary-1 ayer thicknesses are the 
integral momentum equation for <S 2 (x) and the integral energy equation 
for <5 3 (x). The displacement thickness 6j(x) is related to the 
momentum and energy thicknesses through the assumed profile U(y) for 
the boundary layer. 

The governing equations for 6 2 (x) and 6 3 (x) are 


d <$2 (x) 
dx 

d<S 3 ( x ) 

dx 


+ 


(2 + H 12 

V u 





S 3 (x) = C D 


(14) 

(15) 


The shape factor H 12 is a given function of the thicknesses S 2 and 

6 3 and the external velocity gradient . The coefficient Cy is 

the local wall shear stress coefficient and the coefficient Cp is 
the energy dissipation coefficient. 


In the laminar layer, Cy is known from the assumed 
boundary-layer velocity profile and equation (14) can be integrated to 
find <S 2 (x) without solving for the energy thickness. The method of 
Holstein and Bohlen (ch. X of ref. 5) is used to integrate this 
equation. Transition is assumed to occur when the external velocity is 
a maximum; that is, where 


dU < x tr> 

dx 


0 


(16) 


Both equations must be integrated in the turbulent region. 
Truckenbrodt 1 s method (ch. XXII of ref. 5) is used for this purpose. 
The coefficients Cy and Cp are given by empirical functions in 
the turbulent region. 
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The section profile drag coefficient is computed by the method of 
Young and Squire (ch. XXV of ref. 5). In this method, the wake 
thickness at infinity is estimated by the empirical formula 

/ U(x T )\ 3 * 2 

62,00 = M X T ) u — ( 17 ) 

and the section drag coefficient, referred to the chord, is 



The Blade Section Boundary -Layer Module computes the skin friction 
coefficient Cf (fig. 7) for use later in the Propeller Loading 
Module. It provides the drag coefficient for propeller 

performance analysis and for the computation of lifting line (compact 
source) loads. The trai 1 i ng-edge thicknesses 6 1 (x- r ) and 
5 2 ( x i* ) are used in the scaling laws for trai 1 i ng-edge broadband 
noi se. 


Propeller Performance 

Computation of the propeller performance depends on a solution for 
the induced velocity field. This induced field at a blade section is 
shown in figure 8. Propeller performance is predicted by the method of 
Lock (ref. 6). The Blade Section Aerodynamic Module gives tables of 
section lift functions C a (c 1 ,a,M). The Blade Section Boundary-Layer 
Module gives tables of a, M). These tables are used in the 

prediction of propeller performance coefficients Cp and Cj for 
given tip and forward Mach numbers M^ and M of the propeller. 

The helical Mach number at a blade section M^” (fig. 9) is known 
from its components. These components are the axial Mach number of the 
propeller M^ and the rotational Mach number at the section 
^Mt, where” M^ is the tip rotational Mach number Rfi/c^. The 
local induced Mach vector M-j must be found to compute the total 
local Mach vector by 

M h = M. + M (19) 

This vector is represented by its magnitude M and the inflow direction 
angle <j> as shown in figure 9. Differential components of lift dL and 
drag dD on the section are rotated through the inflow angle <p to 
give the differential thrust force dT and torque force dQ. These 
differentials are then integrated over the blade length and converted 
to performance coefficients: the torque or power coefficient Cp 

and the thrust coefficient Cj. The two components of the induced 
Mach vector are supplied by solving the two equations for the change of 
momentum in the far wake of the propeller. The increase in axial 
momentum, the value downstream minus the value upstream, is equal to 
the thrust force. The increase in angular momentum is equal to the 
propeller torque. These balance equations are 
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2M sin <p (M sin <}> - HJ F(£ 1} x) = j M 2 cr( )[C^ cos <}> - C d sin cj>] (20) 


2M sin <f> ( 5 ^ - M cos cf>) F(g ls x) = t)- m2ct ( ) C c £ sin <f> + c d cos <t>3 (21) 


The factor of 2 on the left side of these equations represents the 
fact that velocities in the far wake are twice those at the propeller 
disk. The factor M sin <j> represents the mass flow through the disk. 
The factor (M sin <f> - M„) in the thrust equation is the 
axi al -i nduced Mach number and the factor (^M^ - M cos <j>) is the 
angular-induced Mach number. The Prandtl circulation function F(£j,x) 
(ref. 7) is derived from vorticity theory. It is defined by 


p 

F = — Arccos 
n 



5i) 


(l + X 2 ) 
X 



( 22 ) 


and is an approximate way of representing Goldstein's circulation 
function (ref. 8). This solution procedure was developed in 1930 by 
Lock (ref. 6). 


Loadi ng 

Blade loads are computed by combining results of the aerodynamics, 
boundary layer, and performance modules as shown in figure 10. The 
lift and drag coefficients are three-dimensional tables in terms of 
span position e >1 , angle of attack a, and Mach number M. The surface 
stress coefficients Cp and Cf are four-dimensional tables in 
terms of surface position and ? 2 , angle of attack a, and section 

Mach number M. The performance analysis produces actual values a(£j) 
and MUj) for the propeller operating conditions M ro and M^. 

There may be an additional angl e-of-attack perturbation a'Uj.t) due 
to small nonuniformities in the propeller inflow. 

When the stress and loading coefficient tables are interpolated 
with these functions, the coefficients become actual time-dependent 


val ues 

C £ Ui,t) = C £ [? 1 ,a(c 1 ,t),M(5 1 )] (23) 

C d U!,t) = C d [? 1 ,a(? 1 ,t),M( ?1 )] (24) 

C (?l»?2’ t ) = C [c 1 ,5 2 ,a(5 1 ,t),M(c 1 )] (25) 

r r 

C f (c 1 ,? 2 > t ) = C f [? 1 ,? 2 ,a(C 1 ,t),M(C 1 )] (26) 

where 

aUj,t) = aU ls t) + a'Uj.t) (27) 
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Replacing the parametric arguments a and M by actual arguments 
reduces the lift and drag tables to two dimensions and the stress 
coefficient tables to three dimensions. If steady loadings are 
assumed, the loading tables are reduced by a further dimension. These 
loading tables are passed to the discrete noise prediction module for 
noise prediction. The line loads C Si {E, l ,t) and C<jUi»t) are 
used in compact source theories and the distributed loads 
CpUiJ^’t) and are used in the general noncompact 

source theories. A similar transformation of the boundary-1 ayer 
thicknesses is used to supply data to the Propeller Trai 1 i ng-Edge Noise 
Modul e. 


ACOUSTIC PREDICTION 


Discrete Tone Noise 


Subsonic Noise .- Discrete tone noise is predicted by the method 
of Farassat (ref. 9). Farassat's equation for the noise of a subsonic 
propeller is 


4irp(x,t) 



dS + / 

bl ade 



r 


dS 


(28) 


The terms in equation (28) are illustrated in figure 11. Pressure is 
computed at a particular observer position x and time t. Two 
integrals are evaluated to give the total pressure. The integrals are 
over the surface area of the propeller blade and the integrands are 
evaluated at the time x when the sound is emitted at the surface 
position y(x). The radiation vector r is the difference between the 
observer and source positions 


r = x(t) - y(x) 


(29) 


The normal vector n and surface area dS are given by 




and the velocity vector is 


-► 

V 



x 


+ 

n 


(30) 


(31) 


The loading vector is 

t = {£ (±C f ? 2 + C p n) (32) 

where the positive sign is used on the lower surface between the stag- 
nation point on the lower surface and the trailing edge. The base 
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vector ? 2 is tangent to the surface and in the chordwise direction. 
A unit vector in the radiation direction r is used to define 


M = r • v 
r 

and 

A -> 

£ = r • £ 

Integration of equation (28) i 
solution for retarded time. The ba 
the distance from the source to the 
propagation time: 


( 33 ) 


(34) 

straightforward except for the 
ic retarded-time equation requires 
observer to be compatible with the 


r • r = |x(t) - y ( x) j 2 = (t - t ) 2 (35) 

In the case of a propeller moving along its own axis, the retarded time 
equation can be reduced to 

A g ( t - t ) 2 + 2B 0 ( t - t) + C 0 + Cj cos t = 0 (36) 

Equation (36) has the appearance of a quadratic equation in (t - t) 
except for the coefficient C = C 0 + C x cos x. It can be shown, 
however, that there is a single real solution x < t to this equation 
as long as the motion of the propeller is subsonic. 

Transonic noise .- Supersonic noise is computed by Farassat's 
collapsing-sphere method. The subsonic equation cannot be used on any 
portion of the blade when M r may exceed unity because of the 
jl - M r singularity. In addition, the retarded time equation has 
multiple roots. The collapsing-sphere method is illustrated in 
figure 12. The collapsing sphere intersects the blade surface in a 
curve called the r-curve. Farassat has shown that 

dS _ dr d x / 07 \ 

1 1 - | ” sine 

The integrals for noise are thus evaluated by choosing a set of 
times x and integrating first along the r-curves with fixed source 
time and then over source time. This method was developed by Nystrom 
and Farassat in reference 10. 


Recent advances .- The time derivative in Farassat's acoustic 

equation increases the computation time because at least two integrals 
must be evaluated to numerically compute the derivative. Numerical 
differentiation also introduces some spurious wiggles in the pressure 
signature which appear as increases in the higher harmonics of the 
transformed signal. Recently, Farassat (ref. 11) has taken the 
derivatives inside the integral for both the subsonic and supersonic 
cases. The differentiation under the integral has shown that the noise 
depends on blade surface curvatures. The full implications of this 
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exciting new result are not completely known at this time and are a 
subject of continuing research. The subsonic method used here has the 
derivatives taken inside the integral. The transonic method has the 
derivative left outside the integral. 


Broadband Noise 


Broadband noise is generated by turbulence convected past the 
trailing edge of the airfoil as shown in figure 13. The mean-square 
pressure spectrum is given by Schlinker and Amiet (ref. 12) as 


<p 2 (x,w)> = / 

Trai 1 i ng 
edge 


, k MzC 
l __x 

4 2tto2 


L 2 |£ y U) S( w ) 


ds 

T 


(38) 


where 


0) 

U 

(39) 

= [x 2 + 3 2 (y 2 + z 2 )] 1/2 

(40) 


and 3 is the compressibility factor ( = (1 - The function 

L is the "effective lift" function derived by Amiet (ref. 12) for an 
observer at retarded position (x,y,z). The function £y(w) is an 
empirical correlation length 


2 . 1U C 


where U c is the turbulence convection velocity which is about 0.8U. 

The spectrum function S(ui) is an empirical function for the 
blade surface pressure spectrum at the trailing edge: 


where 


and 


S(u) = (2 x IQ' 5 



2 

u 


F (u) ) 


“ = — 

F(u) = 33 .28u (1 - 5. 49 to + 36. 7w 2 + 0.151w 4 ) -1 


(42) 


(43) 


Near-Field Effects 

Two effects are present in the near field of the propeller which 
may significantly alter the received noise. These effects are the 
scattering by the aircraft wings and fuselage and the refraction by the 
boundary layer on the surface of the aircraft. 
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Scatteri ng .- The scattering effect is illustrated in figure 14 for 
a cylindrical fuselage. The free-field levels on the surface of a 
cylindrical fuselage are calculated by one of the previously described 
methods for discrete noise. This incident field p-j(u),e,x) is 
transformed by 


1 2* » -i(k x x+me) 

M“’ ra ' k x> = 27 J n J e 

U - 00 


P i (w,e,x) dx d 0 


to give the incident field in a wave-number space. The solution for 
the total pressure on the fuselage surface can then be found by 
superimposing the general solutions for incident and scattered 
cylindrical waves such that the boundary condition 


3 p(oj,m,k x ) 

3r 



(45) 


is satisfied. The result for the total surface pressure can be given 
by the transfer function: 



P t (w,m,k x ) = 2T s U,m,k x ) p i (u,m,k x ) 

(46) 

where 

iirk a (l)' -l 



T s (“>r".k x ) = - 2 H m (k r a > 

(47) 


k r ■ C “ 2 - k x 2 ) 172 

(48) 


The factor of 2 in equation (46) is shown explicitly to represent 
the effect of pressure doubling. With this form, the transfer function 
for scattering approaches unity for high frequencies as may be seen 
from the asymptotic forms for the Hankel functions. Following the 
computation of total pressure p t in wave-number space, an inverse 
transform is used to find the surface pressure as a function of 
position ( e,x) . 

Refraction .- The boundary-layer velocity profile alters the sound 
pressure level on the surface of the fuselage by turning the waves 
propagated upstream away from the surface and turning waves propagated 
downstream into the surface. This effect is depicted in figure 15. 

This refraction effect is small for aircraft in low speed flight but 
becomes significant at the higher subsonic Mach numbers. 


A simple model of the refraction effect neglects scattering and 
uses the two-dimensi onal wave equation in a sheared flow as the basis 
for finding a transfer function. The sheared-flow wave equation is 

(» - Uk x ) ^ + 2 ^ k x + („ - Uk x )[(„ - Uk y )2- kkjp = a (49) 
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This equation can be integrated, given the boundary-layer velocity 
profile U(y), from the surface where y = 0 to the edge of the 
boundary layer where y = 6. Initial conditions at the surface are an 
assumed unit pressure and zero pressure gradient: 


p(0) = 1 

(50a) 

dp(°'L = 0 
dy 

( 50b) 


Integrating equation (49) with initial conditions (eqs. (50)) gives the 
pressure p 6 and velocity v 6 at the edge of the boundary layer. 

These results must be scaled to match the pressure in the known 
incident wave field. The process of matching the incident field 
produces a transfer function for refraction: 


T r ( “> k > k x ) = k y P 6 U> k x ) - U - Uk x ) v 6 (a,,k x ) 


(51) 


Numerical integration of the sheared flow equation breaks down at the 
regular singular point <u - Uk x = 0, which corresponds to the 
coincidence of wave speed with flow speed in the boundary layer. 
Special techniques beyond the scope of this paper have been used to 
integrate the sheared-flow wave equation (49) in these cases. 


Combined effects .- The combined effects of scattering and 
refraction may be found by integrating the sheared-flow equation in 
cylindrical coordinates to find the cylindrical wave pressure p 5 
and velocity v 5 at the edge of the boundary layer. Matching the 
external field then gives a combined transfer function: 


k h| 2 ) (z) 


T ( 


sr ' 


j,m,k x ) = 


r r • irZ 

{[ _1 y - 


m 


;nrn 


z )]k r p 5 - [|- 


m 


Z) 


H ( 2 ) (z)' 
m - 


- Mk 


(52) 


where 

z = k p (a + 6) 


(53) 


Far-Field Effects 

Noise propagated to the far field is modified by the effects of 
atmospheric attenuation and ground reflections and attenuation. Since 
phase information tends to be lost after propagation over long 
distances, far-field noise effects are applied to the mean-squared 
pressure spectrum of the noise rather than to the pressure itself 
(ref. 2). 

Atmospheric attenuation .- The effect of atmospheric attenuation is 
computed by the ANSI method described in section 3.1 of reference 2. 

The mean-squared pressure on a source sphere of fixed radius r s is 
reduced by the spherical spreading effect and atmospheric attenuation. 
The transfer function for attenuation T a is a decaying exponential 
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function exp(-2ar). The attenuation rate a is a function of 
frequency. The lower frequencies are dominated by the effect of 
nitrogen relaxation. Mid frequencies are usually dominated by oxygen 
relaxation effects and the higher frequencies are dominated by the 
classical absorption effects of conductivity, viscosity, and rotational 
modes of molecular vibration. The ANSI-proposed standard method is 
used for the calculation of attenuation effects for both standard and 
nonstandard conditions. 

Ground effects.- The far-field noise is reflected and attenuated 
by the ground. This effect may be represented in the transfer function 
T q developed by Pao, Wenzel, and Oncley in reference 13. This 
transfer function is based on the complex reflection coefficient 
Re 1 <f ) for a spherical wave over an impedance plane. The factor R is 
the magnitude of the reflection which is equal to or less than unity. 
The factor <)>• is the phase shift between the reflected and incident 
waves. The magnitude of the noise at the observer depends on the 
difference in the lengths of the direct ray path and the reflected or 
image ray path. The magnitude of the received noise is diminished 
slightly by the loss of coherence of the direct and reflected signals. 
This coherence loss also depends on the path length difference. 


COMPARISONS WITH EXPERIMENTS 
Subsonic Propeller 

Noise predicted by this system has been compared with noise 
measured during the flight of the twin-engined aircraft shown in 
figure 16. The aircraft was fitted with a wing-mounted microphone boom 
which could be moved to measure the noise in front of and behind the 
propeller plane. Noise was also measured on the ground with level 
flights over the measurement point. 

Propeller performance .- An intermediate check of the prediction 
system was made by Block and Martin (ref. 14) by using a 1/4-scale 
model of the twin aircraft propeller on a propeller test stand (fig. 

17) at the Langley Research Center. 

Computed and measured power coefficients are shown in figure 18 as 
a function of advance ratio J = V^/nD and 3/4-span pitch setting 
6 75 . The predicted power coefficient is near the measured values at 
the lower pitch settings but rises above the measured power at the 
higher settings where it is believed that the propeller is in a 
condition of partial stall. Similarly, the predicted thrust 
coefficient Cy (fig. 19) is above the measured value with fair 
agreement at the low pitch settings and poor agreement at the high 
pitch and thrust values. 

Effective pitch.- Pitch was not measured during the flight of the 
aircraft when the noise was measured. In order to find an effective 
value for the pitch setting b >75 to use in noise predictions, the 
measured power in flight was used to compute the power coefficient 
Cp. An effective pitch setting ( B . 75 ) e f f was ttien f° unc * suctl 

that the predicted power coefficient matched the measured value. This 
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effective pitch will be less than the pitch setting on the 1/4-scale 
model for the same Cp. 

Propel 1 er noi se . Near-field noise measured in flight is shown in 
figure 20. The measured noise data were sampled with data blocks of 96 
time points per shaft revolution which gave 32 points per blade for the 
three-blade propeller. Fifty blocks of data were ensemble averaged to 
find the mean signal p(t) shown in figure 20. Signal number 1 is for 
one-third of a propeller revolution. Signals 2 and 3 which complete 
the revolution are similar. The standard deviation a(t) was nearly 
constant, indicating that the measured noise could be decomposed into 
discrete and random parts as follows: 

P (t ) = p(t) + p ' (t ) (54) 

where p'(t) is a stationary random signal. Sound pressure spectra 
were generated for each data block and for the ensemble-averaged data. 

The discrete spectrum |p(w)l 2 and the random spectrum I P 1 ( co ) I 2 are 
added to give the total 1 spectrum |p(w)| 2 . 

PM| 2 = |pM| 2 + |p'(w)| 2 (55) 

The spectral data shown in figure 21 are relative to the overall 
mean-squared pressure 

1 T . 

<p 2 > = lim Yf / P 2 (t) dt (56) 

T+oo 

The overall mean-squared pressure in figure 21 is within 1.8 dB of the 
measured value. This agreement may be due in part to the matching of 
power through the effective pitch setting procedure. 

The predicted spectrum agrees well with the data for the first few 
harmonics. Starting at about the 5th harmonic, the measured discrete 
spectrum shows a cyclic pattern suggestive of cancellation and 
reinforcement which is not predicted. This pattern gives levels from 
the 10th and 15th harmonics which are more than 20 dB above the 
predicted values. A significant part of this underprediction is 
probably due to the omission of unsteady loads in the analysis. 

Despite these high-frequency errors, the predicted A-weighted sound 
pressure level will be in good agreement with the measured values 
because of the rapid decay of the levels with frequency. 

The error in the discrete spectrum may be due to any of several 
effects, including the obvious possibility of an error in the 
prediction method. Scattering from the wing and fuselage was not 
included in these computations. Computations are made for one 
propeller and there may be a small contribution from the second 
propeller which contributes to the error. 
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In making the predictions, it was assumed that all blades are 
identical, whereas, in practice, there are slight bl ade-to-bl ade 
differences which may result in the reinforcement patterns seen here. 
There is regrettably no way to examine these possible sources of error 
without further detailed experimentation. 

Far-field noise from the flyover of the twi n-propel 1 er aircraft is 
shown in figures 22 and 23. The limited data shown in figure 22 are 
for flyovers with the reception angle e = 90° being the directivity 
angle when the sound is received. The prediction is 5 to 6 dBA above 
the average of the measured levels. The reason for the overprediction 
is shown in figure 23. The level for the first harmonic, which 
dominates the predicted source spectrum, is indicated by the data to be 
less than the second harmonic. The reason for this effect is unknown 
and a much larger data set should be examined before any conclusions 
are drawn. Flyover noise data for propellers typically scatter over a 
range of about 10 dB (see ref. 15) so that statistical methods must be 
used to assess the accuracy of predictions. 


Advanced Turboprop 

Predictions have been made for the eight-bladed advanced turboprop 
propeller designated as SR-3. Since this propeller was operating at 
transonic conditions, the elementary aerodynamics and loading methods 
described earlier could not be used. A loading distribution was 
estimated by the propeller manufacturer and this distribution was 
adjusted such that the computed power matched the power measured during 
the flight test of the propeller. 

The propeller was tested in flight atop the JetStar aircraft shown 
in figure 24. Performance data were also obtained in the wind tunnel 
as shown in figure 25. Microphones were mounted in the surface of the 
fuselage along a line under the axis of the propeller. Measured noise 
data were reduced by using ensemble-averaging techniques as described 
earlier. Only the results of the discrete noise data and predictions 
are shown herein. 

Figure 26 shows a comparison of the measured and predicted noise 
at the cruise Mach number of 0.80. The helical Mach number of the 
blade tip is 1.13 for this example. Levels shown are for the blade 
passing harmonic which, remembering the subsonic results, is expected 
to give the best agreement with the measured data. Nystrom and 
Farassat's PROPFAN program (ref. 10) was used to predict the free-field 
noise at the top of the boundary layer. The boundary-layer profile was 
measured with rakes mounted on the fuselage of the JetStar. The 
two-dimensional transfer function for the boundary layer was then used 
to predict the noise on the surface beneath the boundary layer. It is 
apparent that the boundary layer has a significant effect on the 
surface noise at this flight speed. The predictions agree better with 
the data both in front of and behind the plane of the propeller. 
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The causes of the discrepancy between predicted and measured 
levels on the JetStar just behind the propeller plane are being 
investigated at this time. The scattering effect given by equation 
(47) will be evaluated first to see how scattering modifies the 
predicted surface pressures in the absence of a boundary layer. 
Following this the combined effects of scattering and refraction will 
be computed by using the transfer function (eq. (52)). 


CONCLUDING REMARKS 

NASA ANOPP system for predicting propeller noise has been 
described in this section. The objective of this system is to predict 
noise directly from the shape and motion of the propeller blades. This 
objective requi res aerodynamic computations to be made as a basis for 
the acoustic computations. The aerodynamic theories used are simple 
classical methods. 

These methods overpredict the performance of the propellers which 
results in an error in the acoustic predictions. An empirical 
correction procedure of matching predicted and measured power was used 
to compensate for the known error in aerodynamic performance. The 
acoustic predictions give good agreement with the low-frequency 
discrete noise when they are adjusted by the power matching procedure. 
It is believed that improved aerodynamic theories will give good 
low-frequency predictions for subsonic propellers without an empirical 
correction for power. 

The accuracy of the acoustic predictions decreases with increasing 
frequency. Errors of 10 dB or more may occur at or above the 10th 
harmonic of the blade passing frequency in the case of subsonic 
propellers. When the propeller has nearly steady loading, these errors 
do not seriously affect integrated measures of noise such as the 
A-weighted sound pressure level because the A-level is dominated by 
lower harmonics. The prediction methods used are believed to be 
adequate for subsonic propellers operating in a tractor configuration 
where the loads are nearly steady. There is a need for further 
investigation of noise from subsonic propellers operating as pushers 
where the unsteady loads increase the levels of the higher harmonics 
relative to the fundamental. 

Near-field effects of scattering and boundary-1 ayer refraction 
significantly alter the noise levels of the fuselage of the aircraft. 
The boundary-layer refraction effect reduces the surface noise in front 
of the propeller and increases the noise behind the propeller. This 
effect is significant for an aircraft with high subsonic Mach number 
such as the 0.80 Mach number envisioned for advanced turboprop-powered 
aircraft. The two-dimensional model of the boundary layer fails to 
explain a significant discrepancy between predicted and measured data 
just behind the plane of the SR-3 turboprop on the JetStar aircraft. 
Three-dimensional refraction and scattering effects are being studied 
in an attempt to resolve this discrepancy. 
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Far-field effects of atmospheric attenuation and ground 
reflections modify the noise measured during the flyover of a propeller 
aircraft. The largest effect is that of ground reflection which 
introduces the "ground-dip" in the measured flyover noise. Available 
theories and methods are adequate for predicting these effects when the 
ground reflection angle is large. Although measurements from outdoor 
flyover tests have a typically large scatter of data, the flyover noise 
from propellers can usually be predicted within a standard deviation of 
about 4 dB. 
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EXECUTIVE SYSTEM 
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Figure 3.- Propeller blade shape in surface coordinate system. 



Figure 4.- Potential flow mapping for blade section. 



Figure 5.- Propeller blade lift and pressure coefficient functions. 



Figure 6.- Airfoil boundary-1 ayer analysis. 
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Figure 7.- Propeller blade drag and skin friction coefficient functions. 
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Propeller induced flow. 
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Figure 17.- Quarter-scale subsonic propeller in test chamber. 
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Figure 18.- Power coefficient of 1/4-scale subsonic propeller. 










DATA PREDICTION 



Figure 19.- Thrust coefficient of 1/4-scale subsonic propeller. 
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Figure 20.- In-flight near-field subsonic propeller noise signal 
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Figure 21.- In-flight near-field subsonic propeller noise spectrum. 
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Figure 24.- JetStar/propfan flight tests 
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Figure 25.- Propfan performance tests. 
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Figure 26.- JetStar fuselage surface noise. M = 0.80. 



10.2 BLADE SHAPE MODULE 


John W. Rawls, Jr. 
PRC Kentron, Inc. 


INTRODUCTION 

A numerical description of the airfoil surface is the first step 
in rotor noise prediction. Virtually every calculation in the process 
of rotor noise prediction depends on the accuracy of the blade surface 
model. The purpose of the Blade Shape Module is to formulate a func- 
tional representation of the blade surface suitable for aeroacoustic 
and aerodynamic calculations. 

Blade sections are generated by cutting planes normal to the pitch 
axis 02 as shown in figure 1. A sectional view looking outward is 
shown in figure 2. The position and orientation of the section are 
given by the coordinates of the leading edge m a and n 3 g and 
by the twist angle ej which is the angle between the chorcl line and 
the n 1 -ti 2 plane. The shape of the section is assumed to be given 
(fig. 3) in terms of upper surface yy and lower surface y L 
values as a function of distance along the chord x c . 

The first step in the modeling of the airfoil surface is to 
transform the airfoil section data supplied as input from a Cartesian 
coordinate system to an elliptical coordinate system defined by an 
inverse Joukowski transformation (fig. 4). The Joukowski transfor- 
mation, used by Theodorsen (ref. 1) in the early 193U ' s to solve for 
incompressible flow over airfoils of arbitrary shape, offers two 
advantages in modeling the blade surface. The most important advantage 
is in the description of the leading edge. The air flowing over the 
blade surface is disturbed by the leading edge, resulting in the 
emission of sound waves. Thus, the description of the leading edge is 
very important in the noise calculation. In the elliptical coordi- 
nates, 20 to 30 percent of the airfoil shape, starting at the leading 
edge, can be approximated by a single ellipse. As a result, inter- 
polation between data points on such a smooth function gives good 
agreement with the actual airfoil shape. 

A second advantage is that elliptical coordinates allow the blade 
section to be represented analytically by a single monotoni cal ly 
increasing parameter. In the Cartesian coordinate system (fig. 3), the 
airfoil section is a multivalued function of the chord. That is, for a 
given value on the chord line, two values of the surface are given, one 
for the upper surface and a second for the lower surface. The ellip- 
tical coordinates eliminate the problem of having separate analytical 
representations for the upper and lower surfaces. The blade surface is 
represented by the function i n the computational space as 

shown in figure 5. 
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After the section data have been transformed to the elliptical 
coordinates, a cubic spline is fit in the chordwise direction. This 
fit is made in the least-squares sense and a weighting may be applied 
to the input points. With the cubic spline interpolation, the blade 
surface is represented by a two-dimensional computational grid so that 
any point on the surface can be determined by specifying two indepen- 
dent surface parameters, §1 and £ 2 . These parameters give the 
spanwise position n 2 = ?i and the elliptical coordinates 0 = 2t t£ 2 
and 

The final step is to transform the surface points to the pitch 
axis coordinate system. This coordinate system is the reference 
coordinate system for the aeroacoustic calculations. 


A 

AF 

AR 

b 


'Si 


S2 


r 

r £ 

x,y 


x t 


v 




ni,n2»n 3 


SYMBOLS 

blade section area, re R 2 
activity factor 
aspect ratio 

Joukowski transformation parameter 

chord length, re R 

number of output spanwise stations 

number of output chordwise stations 

blade length measured from axis to tip, 
m (ft) 

radius of curvature, re C 

leading-edge radius, re C 

surface coordinates in complex z-plane 

abscissa of leading edge in complex 
z-plane, re C 

abscissa of trailing edge in complex 
z-plane, re C 

blade volume, re R 3 

basis function coefficient 

complex plane of airfoil section represented 
as a near-circle 

pitch axis coordinates, re R 
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,e > 02 ,e > 03 ,e 




0 

0 T 


?2 
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Subscri pts : 

L 

U 


coprdi nates of origin of elliptical axis 
system in pitch axis system, re R 

coordinates of leading edge in pitch axis 
system, re R 

elliptic angular coordinate, 2 tt? 2 » rad 

blade twist angle measured positive clockwise 
looking outboard, rad 

spanwise surface coordinate, 0 <_ gj _< 1 

chordwise surface coordinate, 0 _< ? 2 £ 1 

basic function 

elliptic system radial coordinate 

lower surface 
upper surface 


INPUT 

Description of the airfoil surface begins by defining a sequence 
of cross sections, each of which is defined in a Cartesian coordinate 
system. The leading edge is the origin, with the chord line being the 
x c -axis. The upper and lower surface coordinates are designated with 
respect to the chord line along with a leading-edge radius. This 
section data is given in percent chord length, so that the actual size 
of the blade section is determined by the chord length C. Figure 3 
shows a typical airfoil section in the input Cartesian coordinate 
system . 

Associated with each blade section is a twist angle, defined as 
the angle between the chord line and the rn-n 2 reference plane. The 
reference coordinate system (the pitch axis coordinate system) is a 
right-handed system with axes m» 02 » and 03 * As shown in figure 1, 
the blade is oriented spanwise along the n 2 -axis, with the origin 
through the center of the hub. The description of the blade surface is 
completed by designating the coordinates of the leading edge and the 
blade length R. 

The input consists of parameters which are functions of span and 
the blade section table which is a function of both span and chord. 

The blade section table is input such that the first point is the 
trailing edge, corresponding to x c = 1 and entering upper surface 
points with x c decreasing to x c = 0 or the leading edge. The 
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process is continued over the lower surface with x c increasing 
until again x c = 1. The trailing edge is entered twice due to the 
discontinuity of the slope at the trailing edge. If the leading-edge 
ordinate yu(x c = 0) is not exactly zero, the chord line will be 
rotated to make y y ( 0 ) = 0, and all other ordinates will be adjusted 
accordingly. 

The interpolation of surface points is accomplished by a 
one-dimensional cubic spline. The output controls specify the node 
locations in the chordwise direction through the parameters N^ 2 and 
52 * 


^2 


n 2 

C(n 2 ) 


^( 12 ) 

0 T (n 2 ) 


Spanwise Functions 
number of input sections 
array of spanwise coordinates, re R 
leading-edge abscissa, re R 
leading-edge ordinate, re R 
chord length, re R 
leading-edge radius, re C(n 2 ) 

blade twist angle measured positive clockwise looking 
outboard, rad 


Blade Section Tables 


Ny(n 2 ) 

number 

' of 

input 

poi nts 

on upper surface 



Xy(n 2 ) 

array 

of 

upper 

surface 

abscissa values, 

re 

C(n 2 ) 

YyM 

array 

of 

upper 

surface 

ordinate values, 

re 

C( n 2 ) 

Wy(n 2 ) 

array 

of 

upper 

surface 

weighting constants 


\(n 2 ) 

number 

1 of 

input 

points 

on lower surface 



x L ( n 2 ) 

array 

of 

1 ower 

surface 

abscissa values, 

re 

C(n 2 ) 

y L (n2) 

array 

of 

1 ower 

surface 

ordinate values. 

re 

C(n 2 ) 

w L (n 2 ) 

array 

of 

lower 

surface 

weighting constants 
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Output Controls 
N £ 2 number of chordwise stations 

5 2 array of chordwise stations, 0 _< £ 2 _S 1 


OUTPUT 

The outputs of this module are twi ce-di f ferenti abl e functions of 
the surface coordinates. These functions are the surface position 
vector n-j(?i,S 2 ), the leading-edge position vector ni 
and the elliptic coordinate system data ni eiZ\), b^J, C(Ci), 
and ey^j). All functions have continuous second derivatives defined 
by bicubic splines. This module also outputs the first derivatives of 
the cubic spline basis functions at the node points. 

The aspect ratio, activity factor, and blade volume computations 
are defined in appendix A. 


AF 

activity factor 


AR 

aspect ratio 


Hi 

number of spanwise stations 



number of chordwise stations 


V 

blade volume, re R 3 


5l 

array of spanwise stations, re R, 0 

< 5i < 1 

?2 

array of chordwise stations, re 2 tt , 

0 < ? 2 < 1 


Functions of Span 


n i , Si) 

leading-edge abscissa, re R 


^3 , Si) 

leading-edge ordinate, re R 


bUO 

Joukowski transformation parameter. 

re R 

C(5i) 

chord length, re R 


n l,e (s i> 

elliptic system origin abscissa, re 

R 

''S.e* 5 !* 

elliptic system origin ordinate, re 

R 

0 T(«i) 

blade elliptic axis twist angle, measured positive 
clockwise looking outboard, rad 

AtCj) 

blade section area, re R 2 
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Functions of Span and Chord 


niUi,? 2 ) blade surface abscissa, re R 

n 3 ( ? 1 » ) blade surface ordinate, re R 

'l>(5 1 »5 2 ) blade surface elliptic radial coordinate, rad 


Basis Function Slopes 
4>ij ( Ci-j ) 0 » J = l,2,...,N ?i ) 

i>2j(c 2 i) ( "* ® J - 1,2, . . . ,n^ 2 ) 

METHOD 

The modeling of the airfoil surface follows three basic steps. 

The first step transforms the section data into an elliptic coordinate 
system through the use of the inverse Joukowski transformation. The 
second step fits the blade surface coordinates with cubic splines in 
the chordwise direction. Finally, the third step transforms the 
surface coordinates from the elliptic system to the pitch axis system. 
Each of these steps are discussed in detail. 



Development of Elliptic Coordinate System 
The Joukowski transformati on , defined by 


z = c + ^ (1) 

maps the airfoil section from the complex z-plane to the complex 
5 -plane, where 5 is given by 

c • be^ 1S (2] 

In the 5 -plane the airfoil shape is nearly circular. Recall from 
complex variables that z = re 1 9 is the equation of a circle. If \|> 
is constant in equation ( 2 ), the airfoil shape in the 5 -plane is 
exactly a circle. In general, however, $ is not constant but is 
rather a function of 9. The values of are very small so that 5 
does not deviate very far from a circle. 

To find the representation of an airfoil section in terms of ijj 
and 9, z can be rewritten as 


z = be^ +i 0 + be '*’ 10 


or in terms of hyperbolic functions as 

z = 2b cosh ij i cos 9 + i 2b sinh ip sin e 
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Since z = x + iy in Cartesian coordinates, the coordinates of the 
airfoil section x,y are given by 


x = 2b cosh ip 

y = 2b sinh \p 

Solving for cos 0 and 
following two relations: 


cos e 
sin 0 

sin 0 in terms of x, 


(5a) 

(5b) 

y, and \p gives the 


cos 0 


x 

2b cosh xp 


(6a) 


sin 0 


y 

2b sinh ip 


(6b) 


Using the identity cos 2 0 + sin 2 0=1 together with equations (6) 
yiel ds 

( * ) 2 + [ y 

l 2b cosh \p J V 2b sinh ^ 

For constant ip, equation (7) is the equation on an ellipse with foci 
located at +;2b. Consequently, lines of constant ^ form ellipses in 
the z-plane. Similarly, a relation for cosh ip and sinh \p can be 
found by solving equations (5) in terms of x, y, and 0 to yield 



cosh * ' 2b X cos e (8a > 

sinh * ■ a/sin 9 < 8b > 

Using the hyperbolic identity cosh 2 i|> - sinh 2 4 ) = 1 results in 

x ) 2 _ ( y 

2b cos 0 / (2b sin 0 

Equation (9) is the equation of a hyperbola with foci also located at 
+2b. Lines of constant 0 form an infinite set of hyperbolas in the 
z-pl ane. 

Figure 4 shows an airfoil section plotted in the complex z-plane 
with the elliptic and hyperbolic coordinates superimposed on the 
airfoil. A description of the airfoil surface is sought in terms of 
coordinates ^ and 0 rather than x and y . 



Transformation to Elliptic Coordinates 

In mapping the airfoil to the complex z-plane, the leading and 
trailing edges must be located in the elliptic coordinate system. Note 
that the derivative of the Joukowski transformation 
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^ = 1 - ( 10 ) 
d £ ? 2 

becomes infinite for ?=0 and zero for ?= j^b. By choosing = b, 

the trailing edge becomes a stagnation point. This choice satisfies 
the Kutta-Joukowski condition, which requires the flow field in the 
neighborhood of the trailing edge to be finite even though there is a 
discontinuity in the slope of the surface at the trailing edge. In the 
z-plane, the trailing edge is mapped to the point 

x t = 2b (11) 

Next, the leading edge must be located such that the remaining 
discontinuity is within the airfoil and, thus, will not cause a 
singularity in the flow field. 


Recall from equation ( 7 ) that lines of constant 4) form ellipses 
in the z-plane. The radius of curvature along 4> = Constant is given 
by 

2b (cosh 2 1 p sin 2 e + sinh 2 4) cos 2 e) 3 ^ 2 , 12 
r “ sinh 4< cosh 4/ ' 


With the leading edge located at e = it, the radius of curvature 
becomes 


2b sinh 2 4 
r z ~ cosh ip 


( 13 ) 


Retaining terms of order ip 2 in the Taylor series expansions for 
sinh 4> and cosh 4>, the radius of curvature is approximately given by 

r £ * 2 bip 2 ( 14 ) 


The location of the leading edge in the z-plane can be found by 
substituting y = 0 into equation ( 7 ) and solving for x to yield 

x^ = - 2 b cosh 41 ( 15 ) 

Again by using the Taylor series expansion for cosh ip, equation ( 15 ) 
becomes 

x £ = -2b(l + |i) ( 16 ) 

or in terms of the radius of curvature the location of the leading edge 
is given by 


10.2-8 


-2b 


( 17 ) 


Now the chord length of the airfoil section can be found from 
c * x t - x e = 4b + 


(18) 


Since the chord length is normalized to unity in terms of section input 
data , 


x t - x t ■ 1 


The Joukowski transformation parameter is given by 

r , 


(19) 


K 1 ' i 

b “ 4 "8 


( 20 ) 


Before the airfoil coordinates can be computed in the c-plane, the 
airfoil coordinates must be transformed to the z-plane by the following 
transformati ons : 


x = x„ + x, 
c L 


( 21 ) 


y = 


( 22 ) 


and 


z = x + iy 


(23) 


Solution for Elliptic Coordinates 

Now that the airfoil coordinates have been transformed to the 
z-plane a solution can be found for the coordinates in the c-plane. To 
solve for the elliptic coordinates, variable z is normalized by 2b 
and ? is normalized by b such that 


z = 


_ X + iy 


7F 


( 24 ) 


and 


C = e 


i|)+i 0 


(25) 
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( 26 ) 


Rewriting equation (1) in terms of z and 4 , 
C 2 - 2z c + 1 = 0 
and solving for 5 yields 


? = z ± Jz 2 - 1 

The negative root in equation (27) can be written as 

1 =z-/^7I 

z + /z 2 - 1 


(27) 


(28) 


Consequently, the second solution to equation (26) is the reciprocal of 
the first solution. The elliptical coordinates for the surface points 
z may therefore be found by taking 


ip+ie = In ^ 


(29) 


It is assumed that the function In 4 is uniquely defined by a branch 
cut on the negative real axis, that is, -tt < arg (In 4 ) < 

If the point is an upper surface point, 0 is always positive, so that 


0 = 9 si gn [arg (In 5 ) ] 

^ = <J> si gn [arg (In 5 ) ] 
where sign is defined as 

sign (A) = 

Similarly, for the lower surface, 0 



(A >_ 0) 

(A < 0) 

is always negative so that 


(30) 

(31) 


(32) 


0 = -0 sign [arg (In 4 )] ( 33 ) 

= -ip sign [arg (In 5) ] ( 34 ) 

For 0 to be monotoni cal ly increasing, 2-n must be added to the value 
of 0 on the lower surface. This corresponds to placing the branch 
cut for In 4 on the positive axis such that 0 _< arg (In 4 ) < 2tt. 


Cubic Spline 

Now that the sections are described in elliptic coordinates, the 
airfoil parameters are fit with a cubic spline. (The interpolation 
procedure is described in detail in appendix B.) Introducing new 
variables 


? 1 02 ^" 


(35) 
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to describe the spanwise direction and 

? 2 = 9/ 2tt (36) 

to describe the chordwise direction, a computational grid is formed. 

The N ^2 values of g 2 are input so that 

0 < % 2 < 1 (37) 

As indicated in appendix B, the variables and £ 2 0,0 not have to 

be evenly spaced. 

The cubic spline interpolation procedure allows for two optional 
boundary conditions, one with zero slope or one with zero curvature. 

The zero slope condition gives a rounded edge on the trailing edge, 
whereas the zero curvature condition gives a sharp trailing edge. A 
cubic spline is fit in the chordwise direction £ 2 for in a 

least-squares sense for the case where the number of input points are 
greater than the number of output points. The fit also allows for 
weighting of specific input surface points. If the number of output 
points is equal to or greater than the number of input points, then the 
input grid is used for interpolation. 


Solution for Pitch Axis Coordinates 

The next step is to find the pitch axis coordinates of the grid 
points Uj,^). All lengths in the pitch axis coordinates are 
referred to the blade length R so that the Joukowski parameter is 
scaled by 


bUi) = C(c x ) b (38) 

The pitch axis coordinates are given by 

x(?!, 5 2 ) = 2b(Ci) cosh i|>U lt £ 2 ) cos 2 ttC 2 (39) 

and 

y(5i*€ 2 ) = 2bUj) si nh ^(? ls ? 2 ) sin 2 ttC 2 ( 40 ) 

Let m e j r >2 e» ancl n 3 e denote the location of the origin of 
the elliptic coordinate system in the pitch axis system. Then 

n lje Ui) = n^Ui) - x £ cos e T (c 1 ) (41) 

n 2je = 5, (42) 
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and 


n 3, e ( ? l) = n 3,*^l) + X £ Sin 0 T ( ^l) 


(43) 


where is the location of the leading edge given by 

x ^(?i) = x(5 1s ? 2 =0.5) (44) 

The components of the surface points in the pitch axis system are given 

by 

niUi,? 2 ) = ni.eUi) + [x A (?i»5 2 ) cos e T (c : ) +y £ (c 1 ,£ 2 ) sin e T (t 1 )j (45) 


n 2 (?i,5 2 ) = (46) 


and 

n 3 Ui,? 2 ) = n 3 , e Ui) + [-x^(?i,? 2 ) sin e T (cj) + y 4 Ui,? 2 ) cos 0-,-Ui)] (47) 


Finally, the blade section area is computed by 

1 Sn 3 

A(?j) = / n 1 (? 1 ,S 2 ) -r=- (5i»C 2 ) ds 2 (48) 

0 d ^2 
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APPENDIX A 


ASPECT RATIO, ACTIVITY FACTOR, AND BLADE VOLUME 


The aspect ratio is defined as the blade length squared divided by 
the blade surface area. The equation used in this module is 


AR 



2 


\ l CUj) d?j 

s i,m 


where f, is the minimum value of 5,. 

1 ,m 1 

The activity factor is a measure of the capacity of the blade to 
absorb power. It is defined as 


AF 


100 000 
32 


1 

J CUj) Ci 3 d?! 

^ 1 m 


The blade volume is computed by integrating the blade section area 
as 


V 


1 

/ A Ui) d£i 
'’1 ,m 
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APPENDIX B 


SPLINE INTERPOLATION PROCEDURE 


This appendix 
procedure and then 
approximation (ref. 


detail, which allows for the 
finite-element approach, the 
extended to two dimensions. 


describes the basic cubic spline interpolation 
extends it to a piecewise linear finite-element 
2). The one-dimensional case is described in 


types of boundary conditions. With the 
interpolation procedure is then easily 


One-Dimensional Cubic Spline Interpolation 

One of the simplest ways to interpolate between two data points is 
to construct a polynomial. Assuming that the slope and magnitude of a 
function at two points x-j and x-j + i are known, a cubic 
polynomial can be constructed which is uniquely defined by 


and 


F ( x ) = F. + FJ (x - X.) + 


AF i + l 

A W 


F‘ 2F 1 
r i + l + " i 

Ax i+1 


(x - x i )'■ 


aF 


-2 


i+1 


Ax i + l' 


F Li * F 1 


AX, 


i + 1 


(X - X.) 


(Bl) 


where 


F ( x i ) = F i 


AX-j 


= *1 “ x i_i 


(B2 ) 
(B3) 


AF i 


= F, 


F i-1 


(B4) 


Adding an additional point, two separate cubic polynomials, that is, 
piecewise polynomials, can be constructed which have continuous first 
and second derivatives at xj . Assume f(x) is a cubic polynomial 


on the interval x-j_i £ x < x-j and 
the interval x-j _< x x-j + j. Further 

information is known concerning f(x) 

g(x) is a cubic polynomial 
assume the following 
and g(x): 


f < x i-P = T-i 


(B 5a) 

■ F Li 


(B 5b) 

f(x, ) = g(x i ) = F i 


(B5c) 
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9< x H 1> = F i+1 
9'< x i+l> " F i+1 


(B5d) 

(B5e) 


Differentiating equation (Bl) twice and evaluating both functions f(x) 
and g(x) at x^ yields 


f"tyr ) = 

1 D Ax, 


and 


9 ,,(X i ) 


AF i 

■3 — — - + 2F! + F! , 

Ax . 1 1-1 


AF i + l 

3 11 - 2F 1 - F 

ax . . , i i+1 


i+1 


Setti ng 

f" (Xi ) = g"(x i ) 

and collecting terms of FI yield 


(B6) 


(67) 


(B8) 


Ax 


F! i + 2 (Ax. 


i+1 1 i-1 


i+1 


Ax i ) 


F! 

l 


Ax i 


F 


i+1 


= 3 


Ax "j 


4,1 i+1 




(B9) 


This equation relates the slope at intermediate nodes to the magnitudes 
at the nodes and the slope boundary condition such that first and 
second derivatives are continuous. Equation (69) can be extended to I 
values of x-j , resulting in I - 2 simultaneous equations. The two 
additional equations are provided by the boundary conditions. Two 
options are considered, one with zero slope on the boundary and one 
with zero curvature on the boundary. 

For zero slope, the two additional equations are simply 

F | = 0 (BIO) 

FJ = 0 (Bll) 


For zero curvature, the two additional equations are found by setting 
f"(xi) = 0 in equation (66) and g"(xD = 0 in equation (B7). 

From equation (B7) an expression for zero curvature on the xj 
boundary is given by 
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(B12) 


f ; 4 2 


aF. 


AX- 


and from equation (B6) an equation for zero curvature on the xj 
boundary is given by 


F! 


111 + F* 
2 I 


3 AF 


I 


2 ax 


( B 13 ) 


I 


Having solved for F! , i = 1,2,.. .1, then any value of F(x) can be 
determined by locating 


x i < * < x i + l 


( B 1 4 ) 


and solving equation (Bl). This procedure, which uses a cubic 
polynomial as the interpolating function and requires continuity of 
first and second derivatives at the nodes, is commonly referred to as a 
cubic spline. 


The cubic spline interpolation procedure can be extended to a 
piecewise linear finite-element approach where F(x) is given by 


J 

F(x ) = E 0. <f> • ( x ) (Bl 5) 

j = l J J 


where J is the number of nodes. The function F(x) is defined in 
terms of a set of coefficients Bj and a set of basis functions 
<t> j ( x ) . The basis functions have the following conditions 





(B16) 


and a linear homogeneous boundary condition such that 


♦j (x x ) = ♦j(Xj) = 0 


(j = 1, 2, ..., J) (B17) 


The index i refers to the node and the index j refers to the basis 
function. For J = 3, note in figure 6 that ^ is unity at x x , 
whereas the remaining basis functions <j> 2 and <L are zero. The 
coefficient B], is therefore the magnitude of F[x) at x 1 . 
Similarly, the coefficients Bj are the magnitude of F(x) at the 
corresponding node. Note, also that the nodes x-j do not have to be 
evenly spaced. 
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The set of basis functions <j>j(x) are piecewise continuous 
cubic polynomials over the entire range of x. Consequently an array 
of magnitudes and slopes must be determined for each basis function 
<j>j ( x ) at each node Xj . This results in two matrices, 

♦ ij = (i, j = 1, 2, .... J) (BIB) 

and 

4>!j = ♦ j (x i ) (i, j = 1, 2, .... J) (B19) 

The magnitude matrix <j>-j j is simply the identity matrix defined by 
equation (B16). The slope matrix is computed by solving the 

following equations for the slopes of each basis function 

[ A ki ] K : j] = [B kj ] (i. J, k = 1, 2, .... J) (B20) 


The slope vector solutions from equation (B20) are collected to form 
the slope matrix. Equation (B20) is constructed for the two optional 
boundary conditions by choosing M = 0 and K = 1 for zero slope or 
M = 1 and K = 0 for zero curvature. The k = 1 row of 
equation (B20) is given by 


M ( 2 <j> ! . + <}> ! , •) 
1J i+l, J 


K *ij * 3 M( Vi,j - v /4x m 




( B2 1 ) 


The intermediate rows of equation (B20) are given by 

AX.,, d>! , . + 2(AX., + AX.) d>! . + AX. d> ! , , . 

i+l M-l,j v i+l v 1 n+l,j 


== 3 


AX.. 

AX 




ax, ^ . 

Ax, ^ 6 ij " 6 i-l,j) 


(i = 2, 3, .... I - 1; j = 1, 2, .... J) (B22) 
and finally the k = J row of equation (B20) is given by 

"(♦i-ij + + % ■ 3H(6 ij - (j . i.v !.:. J j) 

(B23) 


With the slope matrix determined, each basis function can be evaluated 
by calculating i such that 

x i < x < x. +1 (B24) 

and evaluating <j>. by 
J 
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-2 


*ij(x - x,) 


(6 i+i,.r 5 ij> 

♦i+u 

AX • , 2 

AX • 

i+l 

1 

5 i + l,j " 6 ij . 



AX. 


i + 1 


AX. 


i + 1 


(x - x. 


(x - x i : 


The coefficients 0j are determined by a least-squares solution 
the equation 


[<0 J (x i r )] {6 d } = {F(xV)} 


j 1 , 2, « • • , J 
1, 2, .... I' > J 


where <t>j( x i') is the matrix of values of the basis function a 
the data points. 


Two-Dimensional Cubic Spline Interpolation 

The one-dimensional interpolation procedure may be extende 
dimensions through the basis functions. The function F(x 1 ,x 2 ’ 
given by 

J i J 2 

F(x 1 ,x 2 ) = l J, $ j i j 2 ^ 1 j ( x i ) ^2 j ( x 2 ) 

j 1=1 j 2 = 1 

The basis functions <J> 1 j(x 1 ) and 4> 2 j ( x 2 ) are piecewise cubic 
polynomials, each a function of only one variable. Therefore, 
basis functions are evaluated in exactly the same manner as in 
one-dimensional case. The coefficients Bj}j 2 are determined 
least-squares solution of a given set of points 
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Figure 1.- Airfoil orientation in pitch axis coordinate system. 



10 . 2-21 


Figure 2.- Typical ai 
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10.3 BLADE SECTION AERODYNAMICS MODULE 



William E. Zorumski 
Langley Research Center 


INTRODUCTION 

Detailed knowledge of the aerodynamic loading on a rotating blade 
is important for the prediction of rotational and broadband noise. The 
Blade Section Aerodynamics Module computes the pressure forces acting 
on the upper and lower surfaces of a two-dimensional airfoil for 
specified angle of attack and Mach number values. This aerodynamic 
output provides input information to the Propeller Loading Module, 
which computes the time dependent aerodynamic loading on the rotating 
bl ade. 

The local flow field on a blade section can be subsonic, 
transonic, or supersonic, depending upon the free-stream velocity and 
location from the hub. Only the incompressible and low subsonic 
computation methodologies are included in this module. 

This module produces a table of the pressure coefficient as a 
function of blade coordinates, angle of attack, and Mach number. 
Theodorsen's method (ref. 1) of successive conformal mapping, which 
yields good comparison with experiments for incompressible flow (i.e., 
M_< 0.4), is used. The Blade Shape Module provides the airfoil shape 
in elliptic coordinates. An application of the inverse Joukowski 
transformation maps this airfoil into a nearly circular curve. An 
inversion of the Theodorsen mapping function transforms this near- 
circle into a perfect circle on another complex plane. Since the 
velocity potential for flow around a circle is well-known, two 
successive transforms produce the flow field around the airfoil. 
Circulation is determined by the Kutta condition at the trailing edge, 
and the leading-edge stagnation point is found by the mapping from the 
perfect circle. The section lift coefficient is found from the 
Kutta-Joukowski theorem. The surface velocity on the airfoil is found 
by the mapping functions and the coefficient of pressure is then found 
by Bernoulli's equation. The incompressible pressure coefficient is 
extended to subsonic flow (i.e., M <_ 0.7) by the simple Karman-Tsien 
compressibility correction. 


SYMBOLS 

A^ ,B n conformal mapping coefficients from ?- to c' -plane 

a radius of perfect circle, re R 

b Joukowski transformation parameter, re R 

C chord length 
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section lift coefficient, — ; 

7 P J 2 S 

oo oo 

section moment coefficient about aerodynamic center, 

positive when nose up, -v- ^ 

\ P v 2 c 2 

C. oo oo 

P - P ro 

local pressure coefficient , — 

i P v 2 

£ 00 00 

complex velocity potential function in ^-plane 


section lift 


Mach number 

section aerodynamic moment 

local aerodynamic pressure 

radial distance in polar coordinates 

blade length, m (ft) 

free-stream velocity 

complex velocity in z-plane 

complex velocity in ^-plane 

physical plane of airfoil in elliptical coordinates 

U,S 2 ) 

blade section angle of attack, measured from abscissa 
of airfoil elliptic coordinate, deg 

circulation, positive clockwise when viewed from hub, 
m-m/s (ft-ft/s) 

= 0 - 2 tt? 2 » rad 

e evaluated at airfoil trailing edge, rad 

complex quantity defining transformation plane 
containing perfect circle c = be ^ o+10 

complex quantity defining transformation plane 
containing an almost circular curve, be ^^ 11 ^ 2 


10 . 3-2 


n 

n 

n 

n 


1 ,c 

3 ,c 

M 

3 ,* 


5 1 


?2 

e 


* 


4 1 o 


p 

Subscripts: 

s 


T 


00 


abscissa of aerodynamic center, re R 

ordinate of aerodynamic center, re R 

leading-edge abscissa, re R 

leading-edge ordinate, re R 

spanwise surface coordinate 

chordwise surface coordinate 

polar angle in transformation plane ?, rad 

blade surface elliptic radial coordinate 

2 TT 

\|) de, rad 
air density, kg/m 3 (slug/ft 3 ) 

stagnation point 
trailing edge 
ambient 



INPUT 

The basic independent variables are the airfoil elliptical 
coordinates ifj(?; 1 ,5 2 ) and arrays of Mach numbers and angles of attack. 


51 

5 2 

b(Ci ) 
C(?i) 


Blade Geometry 

array of spanwise surface coordinates 
array of chordwise surface coordinates 
leading-edge abscissa, re R 
leading-edge ordinate, re R 
Joukowski transformation parameter, re R 
chord length, re R 

blade surface elliptic radial coordinate 
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Output Controls 

N number of angles of attack 

a 

N,. number of Mach numbers 

M 

a array of angles of attack, deg 

M array of Mach numbers 


OUTPUT 

This module produces tables of Cp as functions of surface 
position, angle of attack, and Mach number. This module also computes 
the section aerodynamic lift coefficients C^, moment coefficients 
C m c , the aerodynamic-center locations, and the leading-edge 
stagnation point. 


5 1 

? 2 
a 

M 

C p Ui,5 2 ,a,M) 

C £ (?i,a,M) 

n 3 

? 2 _(?i,a,M) 

9 ^ 


array of spanwise surface coordinates 
array of chordwise surface coordinates 
array of angles of attack, deg 
array of Mach numbers 
local pressure coefficient 
section lift coefficient 

section moment coefficient about aerodynamic center 
abscissa of aerodynamic center, re R 
ordinate of aerodynamic center, re R 
leading-edge stagnation point, rad 


METHOD 

This module computes airfoil section aerodynamics for subsonic 
flow. Theodorsen's method (ref. 1) is adopted for low subsonic flow 
with a pressure compressibility correction for high subsonic flow. 


Conformal Mapping 

The Theodorsen method starts with an airfoil located in an ellip- 
tical coordinate system in the z-pl ane. Under an inverse Joukowski 
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transform, an ellipse in the z-plane of the airfoil will be mapped onto 
a circle. Thus, an airfoil will be mapped onto a s'-plane as a 
near-circle. The Joukowski transformation is given by 

z = c' + b 2 / C * (1) 

A second transform will map the near-circle onto another plane, 

C-plane, as a perfect circle. 

Letting 


_ bg'K ?i» E2) + 2"i S 2 (2) 

be the near-circle and 

c = be** e ie (3) 

be the perfect circle with be^° as its radius, a general trans- 
form from the near-circle to the perfect circle can be written as 

= C exp ( l (A n + iB n )/s n ) (4) 

n=l 

Combining equations (2) and (3) gives 

£* = : ^"^o) e i (2 tt^ 2 - 0 ) (5) 

Equations (4) and (5) yield 

oo 

^ + i (2 -n-C 2 " 0 ) = I ( A n + (6) 

n=l n n 

A general point on the c-plane can also be represented in polar 
coordinates as 


i e 

t, - re 

Hence, equating real and imaginary parts in equation (6) gives 


(7) 


* - + = I 

o 

n=l 


A B 

cos ne + sin ne 
_n 

r r 


( 8 ) 


2ttc 2 - e = l 


n=l L 


B A 

— cos ne sin ne 

„n „n 

r r 


(9) 


10.3-5 


In equation (8), 4> 0 is the zeroth term of Fourier series expansion 

for <J>; thus 

1 f 2ir 

to =h <M 0 (10) 

c J 0 

That is, the radius of the perfect circle is an average of the radii of 
the near-circle. We further have 


A 

r 

B 

r 


n 

n 


£ 

n 


\ 

TT 

I 

ir 



ip cos no de 


sin no do 


(ID 

( 12 ) 


Equations (8) through (12) constitute a system of equations which under 
iteration yields 2 tte; 2 ~ 0. Theodorsen adopted the symbol e to 
denote the shift in arguments going from the to s-plane, that is, 

e = 0 - 2 tt £ 2 (13) 


Iterative Solution for e(c 2 ) 


From the tabulated values of \p and £ 2 » repeated iteration will 
yield new values of ij; and 0, or equivalently, the difference of 
0 - 2 2 • iteration procedure is 

1. Assume an initial approximation, e(s 2 ) = 2 -rr? 2 


l f 

2. Evaluate ip 0 = I iJj(£ 2 (0)) d0; note that numerical 

11 Jo 

integration will involve uneven increments of AO because 
0 is a function of £ 2 which has even increments A e , 2 

A , /'2tt 

3. Evaluate — = - I ip( ^ 2 ( 0) ) cos n0 d0 and 

r Jo 

2ir 

iK? 2 ( 0)) sin n0 d0 



4. Evaluate 2 tt? 2 




cos ne 



sin n0 


5. Loop from step b with 0(5 2 ) = 2 tt£ 2 + e(£ 2 ) 


6. Iteration stops when e(s 2 ) converges 
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Coefficient of Lift 


The flow field about the perfect circle is used to find the 
pressure distribution on the airfoil section. An arbitrary point on 
the surface of the perfect circle is 


l 0 

? = ae 

where 

a = be^° 

The velocity potential in the ?-plane is 


FU) = V ro ( ? e' 1a + y e la ) + In (- e 


C -1 a 


(14) 


(15) 


(16) 


where a is the angle of attack as measured from the real axes 
(airfoil chord) in both the ?-and z-planes. The complex velocity is 
determined from the velocity potential by 


W( c) = — = V (e -1 ° - — e 1 a ) + i— 

d? “ C 2 2 7T^ 


(17) 


The stagnation points are determined by the condition that the complex 
velocity is zero on the perfect circle: 


“ (c) " C |=a 
This is equivalent to 


0 


(18) 


sin ( 0 - a) + 


4iraV 


(19) 


The Kutta condition requires that the complex velocity at the trailing 
edge in the z-jplane be finite. This is equivalent to requiring the 
trailing edge in the ?-plane to be a stagnation point because the 
Joukowski transformation is singular at the trailing edge. At the 
trailing edge, e = so that the circulation r is 

r = AiraV^ sin (a - e^.) (20) 

The leading-edge stagnation point is a solution to equation (19) with 
the circulation given by equation (20). The relationship is 

sin (e - a) + sin (a - e^) = 0 (21) 
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Of the two solutions to the arcsin function, one is the trai li ng-edge 
stagnation point 

e s = <= T (22) 

The leading-edge stagnation point is found from 
sin (e $ - a) + sin (ir - (ct - e-j-)) = 0 

or 

9 s = tt + 2a - (23) 

The section lift is given by the Kutta-Ooukowski law 

4 = PVJ (24) 

so that the coefficient of lift for the section becomes 

= 2t^- S in ( a - e T ) (25) 


Coefficient of Pressure 


The coefficient of pressure is found from Bernoulli's equation and 
the complex velocity in the z-plane. Bernoulli's equation gives 


C 


P 



(26) 


To adjust for compressibility effects, the Karman-Tsien pressure 
correction formula 


C 


P 




incompressible 


(27) 


is applied, and the complex velocity W(z) is found by the transform 


W(z) _ w dz " 1 W( e) 

V ^d?' ; l ds' ] V 
00 00 


(28) 


The independent variable in equation (28) is s', the complex position 
in the Joukowski plane. Both r, and z are given functions of ?' 
as discussed earlier. The derivative of z with respect to ?' is 


dz 

d?' 


U 1 - b)U' + b) 
U ') 2 


(29) 
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so that the expression for the complex velocity becomes 


W(z) .. 


( t - - bn; - .bi ut‘*i -°< b >) 


w(c)/v. 


( 30 ) 


The complex velocity W(c) in equation (30) is 


W(g) - p -io 
V e 


1+2 sin (a - Gj-) 


iae 


i a 


/. i a\ 

2 " 

1 iae 1 


l < y 



and for points on the airfoil surface 


(31) 


dc _ 1 + e' c 

IF ~ 1 - T? F 


(32) 


Primes denote differentiation with respect to £ 2 * For points near 
S* = ±b, equation (30) must be evaluated in the limit as I 
approaches ±b. 


Trail ing edge .- First consider the trai 1 i ng-edge point where 
= b and ? = ?(b) = ^ . Near this point. 


and 


dW(c T ) 

W( c) = W( c T ) + — ( c - c T ) 


(c 


_ \ = dg(b) 
? T ; ~d F 


(I 


b) 


(33) 


(34) 


Since W(^) = 0 by definition, 

(U' - b| « b) (35) 

where 


W(z) 

V 


(c') : 


: d dXl 

b d^ 1 


dg(b) 

d? 1 


2 dW( c T ) 
'L ~ df - 


I dW(g) 

v m d s 


-2ia 


sin (a 



(36) 


Leading edge .- For points on the leading edge near both 
= -b and the "stagnati on point I, an expansion around I = -b 
gi ves 
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dw[c(-b)] 

w(?) = W [c( -b) ] + ^ [c - c(-b) ] 


( 37 ) 


with 


[c - c(-b)] - pp-U' + b) 


(38) 


The expansion for WU) in terms of is then 


WU) = W[c(-b)] 


, dw ^- b) ] drf-b), , 


d? 




^U' + b) 


If the stagnation point is near -b; that is, if 


(39) 


+ b 
s 


« b 


then equation (39) may be used to solve for W[?(-b)] and the 
resulting expression for WU) becomes 


WU) = 


dW[?(-b) ] 


d? 


de(-b) r , 
"de 1 




The expression for W(z) is then 


w(z) _ u 1 ) 2 de(e') de(-b) i_ dW [ g( ~ b) J ’ g s 

v m ? - b de' de' V d? ?' + b 


(40) 



(41) 


The last term in equation (41) is of order unity as long as 
Je^ + bj << b in accordance with the assumption for W U). The field 
point e 1 must also be near -b in order for expansion (39) to be 
used. The two conditions are shown in the following sketch: 
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The largest side of the triangle (-b,e‘ ,?') in the Joukowski plane must 
have a length that is small compared with s b, the approximate radius of 
the near-circle in order for equation (41) to be used. Since two sides 
of the triangle determine the third side, we place the conditions 
on |c‘ + b| and | ^ + b| . The latter condition limits the angle of 

attack a. 


+ 


V i5 2s . 

e + 


« 1 


(42) 


The angle s(c 2 ) is typically small everywhere and iK? 2 ) ls small, 
especially near the leading edge where stagnation is likely to occur. 
With the use of equation (23), condition (42) is therefore approximated 
by 


•e i2 “ + 1 


« 1 


(43) 


or 

|ct| << 1 

Consequently, equation (41) may be used near the leading edge when a 
is small. 

Based on the foregoing discussion, the following sequence of 
computations are used to compute Cp(? 1 ,c 2 ,a). The steps are the 
same for each doublet of values (gj.a) and are as follows: 

1. Compute a by equation (15) 

2. Find the stagnation point e, 2 ,s by equation (23) 
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3. Find ^ by equation (2) 

4. Begin a loop on increasing values of g 2 

5. Compute £* by equation (2) 

6. Compute 0 by equation (13), interpolating on e(£ 2 ) if 

necessary 

7. Compute ? by equation (14) 

8. Compute by equation (32) 

9. If | ?' - b « b, then 

Compute by equation (32) 

Compute y~ — by equation (36) 

Compute by equation (35) 

00 

10. If ?'+bJ < < b, then 

Compute ■ t- by equation (32) 

If ?'+b < < b, then 

s 

Compute y~ by equation (36) 

oo ^ 

Compute by equation (41) 

oo 

Else give warning that is out of range for this 

angle of attack and station 

11. Else compute — by equation (31) 
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12. Compute ^ Z — by equation (30) 

oo 

13. Compute Cp by equations (26) and (27) 

14. Continue looping steps 4 through 13 until all £ 2 values are 

compl eted 

Calculation of C m c and Aerodynamic Center 

The aerodynamic center is defined as the chordwise location where 
dC ro,c = 0. For thin airfoils, its location is observed to be about 
da 

quarter-chord from the leading edge according to the following 
transl ati on: 

r 'l,c = "1,1! * b C0S ®T 

"3,c " "3,1 ' b si " 6 T 

The moment about the aerodynamic center can be computed from the 
pressure coefficient C U 2 ) and the complex coordinate x, by the 
integral p 


/' C p (S 2 ) «e(fb) |f 2 d 52 

where the aerodynamic center is taken at x = -b and Re denotes the 
real part of the complex number. 

REFERENCE 

1. Theodorsen, T.; and Garrick, I. E.: General Potential Theory of 

Arbitrary Wing Sections. NACA Rep. 452, 1933. 
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10.4 BLADE SECTION BOUNDARY-LAYER MODULE 


Donald S. Weir 
PRC Kentron, Inc. 


INTRODUCTION 

The Blade Section Boundary-Layer Module computes the two- 
dimensional boundary layer on airfoil sections using integral 
formulations for the boundary-1 ayer thicknesses. The methods here are 
taken directly from Schl ichting' s Boundary-Layer Theory (ref. 1). The 
blade surface arc length is derived from data supplied by the Blade 
Shape Module and the external velocity distribution is supplied by the 
Blade Section Aerodynamics Module for each section angle of attack and 
Mach number. The method of Holstein and Bohlen is used to integrate 
the momentum equation in the laminar regime. Transition is assumed to 
occur at either x m , where the velocity is maximum, or at a trip 
location x^ r supplied by the user, whichever is the lesser. 
Truckenbrodt 1 s method is used to integrate the energy equation from the 
transition point to the trailing edge. Error messages are given if 
separation occurs in the laminar regime. The profile drag is computed 
by the method of Squire and Young. The module outputs the displacement 
and momentum thicknesses, the skin friction coefficient, and the 
profile drag coefficient for the section. 


b 

C 

C D 


Cf 

Ct 

Coo 

922 

H 


SYMBOLS 

Joukowski transformation parameter, re R 
chord length, re R 

turbulent energy dissipation coefficient. 


00 

ambient speed of sound, m/s (ft/s) 

metric coefficient for surface arc length, re R 

modified shape factor 


l * K£r) 


local skin friction coefficient. 


local skin friction coefficient. 
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M 

R 

R 2 

r 3 

Roo 

U 

u 

X 

a 

a(H) 

B(H) 

5 

6 2 

V 

51 

5 2 

P 

T 

* 


Mach number of blade section 

blade length measured from axis to tip, m (ft) 

U6 2 

Reynolds number based on momentum thickness, — — 

U <5 3 

Reynolds number based on energy thickness, — — 

cj* 

Reynolds number, 

potential flow velocity, re c^, 

velocity inside boundary layer, re c ro 

arc length measured from stagnation point, re C 

angle of attack, rad 

skin friction function 

energy dissipation function 

boundary-layer thickness, re C 

6 

displacement thickness, / (1 - tt) dy, re C 

momentum thickness, / (1 ~ jj) dy, re C 

6 

energy thickness, / [1 -(tt) 2 ] tt dy, re C 

0 u u 

kinematic viscosity 
spanwise surface coordinate 
chordwise surface coordinate 
density, kg/m 3 (slugs/ft 3 ) 
skin friction stress, N/m 2 (lb/ft 2 ) 
blade thickness function, rad 


10.4-2 


Subscripts: 


l lower surface 

m maximum velocity point 

s stagnation point 

t trailing edge 

tr trip location 

u upper surface 

°° ambient 


R 


oo 


b(?i) 

cu x ) 

? 2,tr,u^^ 
? 2,tr, SL^ 1 ^ 


?2 

<J>Ui,S 2 ) 


INPUT 

User Parameter 

Reynolds number based on blade length, speed of sound, 
and ambient viscosity, Rc m /v 

Functions of Span 

array of spanwise surface coordinates 
Joukowski transformation parameter, re R 
airfoil section chord, re R 
upper surface trip location, rad 
lower surface trip location, rad 

Functions of Span and Chord 
array of spanwise surface coordinates 
array of chordwise surface coordinates 
blade surface elliptic system radial coordinates 


Functions of Position, Angle of Attack, and Mach Number 
a array of angles of attack 

M array of Mach numbers 
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5 c (?i»a>M) stagnation point location, rad 

^ 9 S 

Cp(c 1 , 52 » a > |Vl ) coefficient of pressure 


OUTPUT 


a 

M 

Cd ( 5 i > ) 
«i f u,t(Ei,M) 

5 2,u,t(Si,a,M) 

6 l,£,t(5i»a»M) 

(S 2,£,t(5i»a,M) 

Cf(5i,?2> a > M ) 


array of spanwise surface coordinates 
array of chordwise surface coordinates 
array of angles of attack, deg 
array of Mach numbers 
profile drag coefficient 

upper surface trai 1 i ng-edge displacement thickness, 
re CUj) 

upper surface trai 1 i ng-edge momentum thickness, 
re CUi) 

lower surface trai 1 i ng-edge displacement thickness, 
re CUi) 

lower surface trai 1 i ng-edge momentum thickness, 
re 0 ( 5 ^ 

skin friction coefficient 


METHOD 

Boundary-Layer Equations 

The integral equations (ref. 1) for boundary layers are based on 
several thickness functions, shape factors, and coefficients defined 
here for later use in analysis and computation. The basic thickness 
functions are the displacement thickness 

6 

«! = / (1 - u/U) dy (1) 

0 

the momentum thickness 
6 

S 2 = f (1 - u/U) (u/U) dy (2) 

0 


C-* 
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and the energy thickness 


S 3 = / (1 - (u/U) 2 ) (u/U) dy 

0 

Reynolds numbers based on these thicknesses are defined by 


K i ” v 


(i = 1, 2, 3) 


Shape factors which give the ratios of these thicknesses are 
6 R • 

H ij = 67 = IT (i.j = !» 3; 1 * J) ( 5 < 

The skin friction coefficient based on local velocity at the edge of 
the boundary layer is 


~ \ v / 

T puV 

oo 

and the skin friction coefficient based on the blade section velocity 


Cf = i 

f i pH 2 c 2 
2 00 


so that 


C f " 2 W C T 


The energy dissipation coefficient is 


VpU 2 c 2 
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The basic governing equations for momentum and energy in the 
boundary layer are 


ds 2 

dx 


+ 


6 2 

(2 + 



and 


d6 3 63 dU 

d>T + 3 U~ dx 2C D 


( 10 ) 


(ID 


The displacement thickness is found from these thicknesses and an 
assumed form of the boundary-1 ayer profile. This profile is different 
in the laminar and turbulent cases and is defined later. 

Arc Length .- The boundary-layer analysis is carried out on both 
the upper surface and lower surface of the airfoil. In this analysis, 
these two surfaces are defined in terms of two continuous arcs measured 
from the stagnation point. The arc length on the airfoil surface is 


dx = ±/ g 22 d? 2 

where the metric coefficient for arc length is 

\2 


( 12 ) 


922 = 4 (£) 1 + OHf) 


(sinh 2 ^ + sin 2 2 tt£ 2 ) 


( £ = Constant) (13) 

The arc length begins with x = 0 at the stagnation point £ 2S . 
Equation (12) is integrated on the upper surface with the negative sign 
in equation (12) being used because the arc length increases as £ 2 
decreases to 0 at the trailing edge. On the lower surface, the 
positive sign of equation (12) is used with £ 2 increasing to 1 at the 
trailing edge. 


of 


Local 

pressure 


velocity.- The local 

“C by 
P 


velocity is found from the coefficient 


u(x) = m n 


(14) 


This velocity is actually a Mach number since the speed of sound c«, 
is used for a reference. The symbol M is saved for the local section 
Mach number which does not depend on chord position. 

Maximum velocity .- Both U(x) and U'(x) are used in the 
computation. After finding U'(x), the maximum point for U is found 
from 
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(15) 


dU(x ) 
v nr 

dx 


0 


(x > 
v m 


0 ) 


Local Reynolds number .- With U(x) given relative to c* and 
6(x) given relative to chord length C, the kinematic viscosity must 
be relative to Cc*,. Thus the local kinematic viscosity is given in 
terms of the input Reynolds number parameter R*, = Rc OT /v by 

v = tttt (16) 


The local Reynolds number is 

R = MCR = — (17) 

oo v 


Laminar Boundary Layer 

The laminar boundary-layer equations are integrated by the method 
of Holstein and Bohlen as described in chapter X of reference 1. The 
integration is carried out with a shape function K(x) defined by 


K(x) 



dIJ 

37 


and a function of the momentum thickness 


(18) 


«2(x) 

z(x) = — — (19) 

Four auxiliary functions of K are used in the integration. These are 
given in table 10.2 of reference 1 and here as table I and figure 1. 

The auxiliary functions are Pohlhausen's shape factor A 

-12! < A(K) < 12 (20) 

the shape factor H 12 

H 12 = fiOO (21) 

the wall stress factor f 2 (K) which gives Cj by 

R 2 C t = f 2 (K) (22) 

and a function F(K) in the following differential equation for Z 

“ = F(K) + K (23) 
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As equation (23) is integrated, K is calculated from 


K ■ z 37 (24 > 

The limit A = -12 indicates separation and the limit A = 12 
constrains the boundary-layer profile such that max [u(y)/Ul < 1. 

Stagnation point .- The value of Z at the stagnation point is 

z <°) ' TO) (25) 


where 


K ( 0) = K 0 = 0.0770 (26) 

The boundary-layer momentum thickness at the stagnation point is 


S 2 (0) = /v7(0T" (27) 

and the displacement thickness is 

6i(0) = f i (K 0 ) 6 2 (0) (28) 


The skin friction coefficient at the stagnation point is 
C f (0) = 0 


(29) 


Laminar range .- Equation (23) is integrated over the laminar range 
until either (1) the trailing edge is reached or (2) A is out of 
range, | A J >12 . 0 . The end of this region is denoted X[_ and saved 
for future reference. In the laminar range 0 < x < X|_, the 
boundary-layer parameters are 


6 2 (x) = / ( vZ ( x ) ) 



(30) 

M x ) = f i(K(x))6 2 (x) 



(31) 

and 




2U(x) v f 2 ( K(x ) ) 

C f (x) = 

M 2 6 2 (x) 



(32) 

The skin friction coefficient is taken to 
surface where 0 < £ 2 < ? 2>s and positive 
5 2 s < ?2 < 1* TFis convention allows the 
the direction of £ 2 everywhere. 

be negative on the upper 
on the lower surface where 
friction stress to be in 


Transition 


The transition from laminar to turbulent flow is influenced by a 
number of factors, but the strongest influence (Ch. XVII of ref. 1) is 
the pressure gradient in the boundary layer. When this gradient is 
positive, the boundary layer is unstable and transition to turbulent 
flow soon follows. It is assumed here that the transition point is the 
point of minimum pressure or the point of maximum velocity. Thus, the 
point of transition x^ r is found by solving 


dU(x 

3x 



(x > 0) 


(33) 


for the transition point x^ r . Since equation (33) is identical to 
equation (15), we use x^ r = x m . 


Turbulent Boundary Layer 

Modified shape factor .- An additional shape factor introduced by 
Truckenbrodt (Cn. XXII of ref. 1) is used in the calculations for 
turbulent boundary layers. Truckenbrodt 's modified shape factor is 
defined by 



(34) 


The velocity profile assumed for the turbulent boundary layer is such 
that 


H 32 

” 12 2 (3h 32 - 4) 


(35) 


so that equation (34) can be integrated. The value of H 32 for a flat 
plate is denoted H 32 » and is the lower limit of the integral in 
equation (34). 


Truckenbrodt 's modified shape factor is then 


/ ^23 

H = G - H - 

where 

Gi = 0.5049 G = 0.5442 

1 oo 


(36) 


(37) 
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are empirically determined constants. Equation (36) may be inverted by 
solving a cubic equation in H 23 . The solution is given in two steps 


<l> - -j Arctan 


Energy method.- Truckenbrodt's energy method is based on 
integrating differential equations for the energy Reynolds number R 3 
and the modified shape factor H. The differential equations are a 
transformation of equations (10) and (11). 


33T ° " 2 (ij (e) r 3 + u *3<f<3. H ) 


The functions $(R 3 ,H) and i'(R 3 ,H) are empirical functions which are 
described later. 

Initial conditions .- The initial conditions for these differential 
equations are found from H 12 = f 3 (K) - AH 12 at the end of the laminar 
range where x = x m = x^p. The correction aH 12 is an empirical 
function relating the laminar and turbulent shape factors and is given 
in table II and figure 2. From equation (35), 


H 32( x tr ) = 3H77 


and, from equation (36), 


H(x tr } = H 32 
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The boundary-layer thickness is 


<5 3^ x tr^ ~ ^ 32 ^ x tr^ <S 2^ x tr^ 


and the Reynolds number is 


Mv' * 


h(\ r 1 u(x tr ) 


( 44 ) 


(45) 


Equation (45) gives the initial condition for R 3 (x) and equation (43) 
gives the initial condition for H(x). 

Empirical functions .- The empirical functions give the effects of 
the skin friction and energy dissipation coefficients in equations (10) 
and (11). These functions are defined in reference 1 as 


*3 



and 


(46) 


*3 


H 

v 


- H 32 ^ J 

H 12 - 1 


(47) 


The shape factor ri 32 is found from equation (39) and H 12 is found 
from equation (35). The coefficients Cy and Cg are defined as 


C 


T 



(48) 


C D 

_ 6(H) 

r 3 


(49) 

a. 

b, and 3 

given as 


a = 

0.268 


(50) 

a = 

0.245(1 - 

2.007 log H 12 ) 1 - 705 

(51) 

b = 

0.2317H 32 

- 0.2644 - (0.896 x 10 5 ) (2 - H 32 ) 20 

(52) 

3 = 

[0.00481 + 

0.0822(H 32 - 1.5) 4 - 81 ]H 32 b 

(53) 


Turbulent separation .- The point of separation for the turbulent 
boundary layer is a function of the modified shape factor H. This 
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factor is near unity at the transition point and generally decreases 
with distance along the airfoil. The condition for attachment of the 
boundary layer is thus 

H > H $ = 0.723 (54) 

The boundary layer is prone to separation at slightly larger values of 

He, referred to as H' where 
b s 

= 0.761 (55) 

The separation point should be set by the user with a parameter 
0.723 < H s < 0.761. 


Profile Drag 

The profile drag or total drag of an airfoil is the sum of the 
skin friction drag and the pressure drag or form drag (ch. XXV of 
ref. 1). The method of Squire and Young presented in reference 1 is 
based on the wake momentum thickness at the trailing edge: 

5 2 i = <s 2 ,u( x u,t) + <5 2,«.( x £,t) (56) 


The profile drag coefficient is given by the wake momentum thickness 
far from the trailing edge as 


C 


d 



(57) 


The method of Squire and Young relates the far wake momentum thickness 
to the trailing edge thickness by the empirical formula 


6 2 



6 2 



(58) 


where U x is the potential velocity at the trailing edge. Since there 
may be small numerical differences between U(xt) on the upper and 
lower surfaces, U x is defined here as the average: 


Ui 


U (x . ) 
u v u ,t 


U i< X £.t> 


2 


(59) 
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Computation Method 



For each section, angle of attack, and Mach number. 


Compute the local viscosity v 
Compute U(s 2 ) from equation (14) 


Compute at the stagnation point ? 2 ,s 

2. Compute the arc lengths on the upper and lower surfaces 

3. Locate the point of maximum velocity x m on the upper and lower 

surfaces. Repeat steps 4 to 7 for each section 


4. Find the initial conditions at the stagnation point x = 0 for the 

laminar functions K and Z 

5. Integrate the laminar equation to the transition point xtr = x m 

or a point x^ r < x m specified by the user, stop the 
integration and issue an error if laminar separation occurs 

6. Find the initial conditions for the turbulent equations at x^ r 

7. Integrate the turbulent equations to the trailing edge or until 

separation occurs as indicated by the limiting shape factor 

H s 

8. Compute the profile drag coefficient. 


REFERENCE 

1. Schlichting, Hermann (J. Kestin, transl.): Boundary-Layer Theory, 

Seventh ed. McGraw-Hill Book Co., c.1979. 
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TABLE I.- AUXILIARY FUNCTIONS FOR LAMINAR 
BOUNDARY-LAYER COMPUTATION 


K 

A(K) 

F(K) 

f i(K) 

f 2 (K) 

-0.1567 

-12 

1.7241 

3.500 

0 

-0.1474 

-11 

1.6257 

3.383 

0.019 

-0.1369 

-10 

1.5229 

3.276 

0.039 

-0.1254 

-9 

1.4167 

3.176 

0.059 

-0.1130 

-8 

1.3080 

3.085 

0.079 

-0.0999 

-7 

1.1981 

2.999 

0.100 

-0.0862 

-6 

1.0877 

2.921 

0.120 

-0.0720 

-5 

0.9780 

2.847 

0.140 

-0.0575 

-4 

0.8698 

2.779 

0.160 

-0.0429 

-3 

0.7640 

2.716 

0.179 

-0.0284 

-2 

0.6609 

2.647 

0.199 

-0.0140 

-1 

0.5633 

2.604 

0.217 

0 

0 

0.4698 

2.554 

0.235 

0.0135 

1 

0.3820 

2.508 

0.252 

0.0264 

2 

0.3004 

2.466 

0.268 

0.0385 

3 

0.2255 

2.427 

0.283 

0.0497 

4 

0.1579 

2.392 

0.297 

0.0599 

5 

0.0979 

2.361 

0.310 

0.0689 

6 

0.0459 

2.333 

0.321 

0.0706 

6.2 

0.0363 

2 .328 

0.324 

0.0721 

6.4 

0.0274 

2.323 

0.326 

0.0737 

6.6 

0.0186 

2.318 

0.328 

0.0752 

6.8 

0.0102 

2.314 

0.330 

0.0767 

7 

0.0021 

2.309 

0.331 

0.0770 

7.052 

0 

2.308 

0.332 

0.0781 

7.2 

-0.0051 

2.305 

0.333 

0.0794 

7.4 

-0.0132 

2.301 

0.335 

0.0807 

7.6 

-0.0203 

2.297 

0.337 

0.0819 

7.8 

-0.0271 

2.293 

0.338 

0.0831 

8 

-0.0335 

2.289 

0.340 

0.0882 

9 

-0.0608 

2.273 

0.347 

0.0919 

10 

-0.0800 

2.260 

0.351 

0.0941 

11 

-0.0912 

2.253 

0.355 

0.0948 

12 

-0.0948 

2.250 

0.356 
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TABLE II.- 


VARIATION OF SHAPE FACTOR H 
AT TRANSITION 


R 2,tr 

AH 12 

1 x io 3 

1.16 

2 x 10 3 

1.22 

5 x 10 3 

1.27 

1 x 10 4 

1.31 

2 x 10 4 

1.34 

5 x 10 4 

1.35 

1 x 10 5 

j 

1.36 


f 
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Figure 1.- Auxiliary functions for laminar boundary-! ayer computation. 
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Variation of shape factor H-.o at transition. 


10.5 PROPELLER PERFORMANCE MODULE 


Ollie J. Rose 
PRC Kentron, Inc. 


INTRODUCTION 

A knowledge of the effective velocity field in which a propeller 
operates is essential in order to predict the blade loading and noise. 

A prediction of thrust and efficiency is needed for flight path 
calculations and for comparison with experiments. The unknown quantity 
required for these purposes is the induced velocity field. Thus, the 
purpose of the Propeller Performance Module is to compute the induced 
velocity field and then the thrust coefficient, power coefficient, and 
efficiency for a given propeller under specified operating conditions. 
The induced flow calculation is performed with the propeller axis 
operating at zero angle of attack. Then, the kinematic effect of 
propeller angle of attack is included by a rotational transformation of 
the total velocity field. 


a i 

a 2 


c 

c 


c 


00 


F 

M 


SYMBOLS 

induced axial velocity component, re ARft 
induced tangential velocity component, re rRft 
blade section chord, re R 
section drag coefficient 

section lift coefficient 

power coefficient 

torque coefficient, = C p 

section torque coefficient 

thrust coefficient 

section thrust coefficient 

speed of sound, m/s (ft/s) 

tip relief factor 
Mach number 
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N 

P 

Q 

R 


r 

T 

V 

V. 


x.y.z 


a 


a 


P 


e 


n 

0 

0 

x 

I 


a 

4 > 

<l> 

ft 


blade tip Mach number, — ^ 

oo 

inflow Mach number normal to disk plane 
number of propeller blades 
power, re p R 5 n 3 

oo 

torque, re pjl 5 ^ 2 
radius of propeller disk, m (ft) 
radial position on blade, re R 
thrust, re p R 4 ft 2 

oo 

effective velocity relative to blade section, re Rft 

resultant velocity of fluid in thrust direction 

resultant velocity of fluid in disk plane in direction 
of rotation 

right-hand reference coordinate system, z in thrust 
di rection 

blade section angle of attack, rad 

propeller angle of attack, rad 

relative error 

propeller efficiency 

blade twist angle, rad 

root pitch setting, rad 

local advance ratio 

average advance ratio 

ambient density, kg/m 3 (slug/ft 3 ) 

blade element solidity, - — 

inflow angle, rad 

blade rotation angle, rad 

angular velocity of blade, rad/s 
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Subscripts: 


partial derivative with respect to M 
partial derivative with respect to a 


9 


a 


INPUT 

The quantities required by this module are grouped into four basic 
classes: blade geometry, flow field, operating parameters, and 

aerodynamic characteristi cs. A uniform or radially varying flow field 
can be specified. 


C(r) 

e(r) 


r 

M z (r) 


N 

M. 


(-) 

T'tii 


Blade Geometry 
number of spanwise stations 
array of spanwise stations, re R 
blade chord array, re R 
section twist angle, rad 

Flow Field 

number of spanwise stations 

array of spanwise stations, re R 

inflow Mach number perpendi cul ar to propeller 
rotation disk 

Operating Parameters 
number of blades 

R 0 

blade tip Mach number, — 

oo 

root pitch setting, rad 

ratio of solidity to tip loss factor at the blade 
ti p 
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a p 

propeller angle of attack, rad 

e 

error parameter 


Aerodynamic Characteristics 

N r 

number of spanwise stations 

r 

array of spanwise stations 

N 

a 

number of angl e-of-attack values in aerodynamic 
tabl e 

a 

array of angl e-of-attack values, deg 

M 

number of Mach number values in aerodynamic table 

M 

array of Mach numbers 

C d (r,a,M) 

section drag coefficient 

C^(r ,a,M) 

section lift coefficient 


OUTPUT 

c p 

power coefficient 

C T 

thrust coefficient 

7 

advance ratio 

n 

propeller efficiency 


local angle-of-attack, rad 

M(r,i|)) 

local Mach number 

V z (r,i|») 

resultant velocity of fluid in the thrust direction, 
re RQ 

V r **> 

resultant velocity of fluid in disk plane in direction 
of rotation, re Rfi 

a 1 (r,^)x(r) 

induced axial velocity, re Rn 

a 2 (r,^) 

induced angular velocity, re n 


inflow angle, rad 
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METHOD 


Induced Flow 

Components of the incoming flow field are provided relative to the 
propeller disk. Radial velocity components and other three-dimensional 
effects are neglected; thus the use of blade element-momentum theory 
together with two-dimensional aerodynamic characteristics is allowed. 

In addition, the principle of momentum conservation is applied to 
annular control volumes and any interference effects between adjacent 
annuli are neglected. A further simplification is achieved by ignoring 
interference effects between successive blades; thus, the theory is 
restricted to axially symmetric inflows and induced velocities. The 
effects of propeller angle of attack are incorporated by a simple 
rotation of the resulting flow field. 

An examination of figures 1 and 2 reveals that the relative 
section velocity at radial station r is given by 

V(r) = {[x(l + ai )] 2 + [r(l - a 2 )] 2 | 1/2 (1) 


where 


A(r) = M z (r)/M t 


The local 

inflow angle 

and angle of attack are given by 




A(1 + a j) 



tan = 

r(l - a 2 ) 

(2) 

and 





+ 

O 

CD 

II 

e(r) - 4 >(r) 

(3) 


The quantities 0 , e(r), and x(r), are presumed known, whereas a x 

and a 2 are to be calculated. 


Equations for a! and a 2 

Consider an annular control volume at station r and of width dr 
as indicated in figure 3. The element of thrust on this annulus is 
given by 

dT = -^-V 2 N C(r)[C cos 4 > - sin <j>] dr (4) 

Similarly, the element of torque on the same annular region is 

dQ = jV 2 N C(r)[C^ sin <f> + C d cos 4 >]r dr (5) 
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On the other hand, linear and angular momentum increase through the 
annulus give the following expressions for the element thrust and 
torque: 


dT = 4irrx 2 a 

i(l + 

a : )F 

dr 




(6) 

dQ = 4ur 3 xa 

2 (1 + 

a 2 )F 

dr 




(7) 

Equating the two expressions 

for 

thrust and torque 

gi ves 



4TrrX 2 a 1 (l + 

a i)F 

1 

2 

V 2 N C(r) 

[C ^ cos 

<f) - sin 

♦ ] 

(8) 

4TTr 2 xa 2 (l + 

ajF 

1 

2 

V 2 N C(r) 

1°Z s1n 

p + cos 

♦ ] 

(9) 


Division of equation (9) by equation (8) gives 

ra 2 C £ sin <J> + C d cos <J> 

Aa i ~ C co~<f> - C d sin 4 > 


Using the notation 


and 


Cq = sin ^ + C d C0S ^ 


C = C cos p - C , sin <)> 


't — r ~ d 

allows equation (10) to be written as 

a, 

Similarly, equation (8) becomes 
4 X 2 a i(l + a i)F = V 2 oC t 


( 11 ) 


( 12 ) 


(13) 


(14) 


The solidity a(r) is defined as 


cr(r) 


N C(r) 
2irr 


(15) 


The tip relief factor F is taken to be unity for the blade 
element-momentum theory. A second option exists to account for 
three-dimensional effects. The Prandtl theory (refs. 1 and 2) gives 
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an approximation for the radial flow effects based on the potential 
flow around the edge of a set of semi -i nf i nite vortex sheets. The 
resulting value of the tip relief factor is 


F = — Arccos S exp 


•N( 1 - r) ( 1 + A 2 ) 1/2 

2x 


Note that F approaches zero with infinite slope as r approaches 
unity. For a nontrivial solution to equation (14), the ratio of a / F 
must remain finite at the tip. If the chord C(r) approaches zero 
with finite slope at the tip, then by L'Hopital's rule, the ratio ct/F 
approaches zero at the tip. The user may specify a nonzero limit for 
the ratio (ct/F)^. . 


Solution Procedure 

Solutions to equations (13) and (14) correspond to zeros of the 
following functions 


Hi(a i,a 2 ) = C Aa j - rC.j.a 2 


H 2 (a 1 ,a 2 ) e 4 A 2 a : (1 + ajF - V 2 a C t (18) 

This system of equations is solved iteratively using a general form of 
Newton's method. The procedure is written in vector form as 


a (n+D 

a i 

a (n+l) 


i(n) ,(n) 
J 1 1 d 12 

3 (n ) ,(n) 

O 2 1 o 2 2 


The elements J-jj, i , j = 1, 2, of the Jacobian matrix are determined 
from equations (17) and (18) as 

8H j A 2 a 2 

J ll(a l»a 2 ) = 3 ^ = C q A ' 1 T~t C0S * (sin ^ C i,a + C0S * c d, a ‘ C t } 

- M sin <p (sin $ C £>M + cos * C d M ) ] 
rAa 2 

+ T 7 — [cos <p (cos d> C„ - sin <t> C . + C ) 

V 1 T T £,a T d,a q 


- M sin * (cos <j> - sin <(, C d>M ) ] 
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^ 12 ( a 1 > a 2 ) = 


n 


3H j Ad ^ 

Si! = r (- C t - —I 51 '" * < Sln + C e ,„ + cos * C d,„ - C t> 

+ M cos <j> (sin <p C M + cos <}> C d M ) ] 


+ n (sin <)> (cos <k C „ - sin * C, + C ) 

V L r l, a r a, a q 

+ M cos (#> (cos <f> C £jM - sin <j> C d M )]} ( 21 ) 

9H 2 

^ 2 1 ( a l » 3 2 ) = "gg = ^ A 2 F ( 1 + 2a : ) - 2V a\ C^. sin <j> 

+ VaA(cos d> (cos d> C - sin d> C , + C ) 

L Y £, a y d, a q' 

- M sin <j> (cos * C M - sin 4 , C d M ) ] (22) 

9H 2 

J 2 2 ( a 1 » a 2 ) = = Vorfsin <i> (cos 4 C - sin 4, C . + C ) 

9a 2 Z,a a »a q 

+ M cos 4 > (cos <f> M - sin 4 > C d iV) + 2C t cos $)] (23) 

With the initial values a^ = a 2 °^ = 0, the iteration defined by 

equation (18) is performed until the following error criterion is met: 

|(a<" +1) ,a<" +1 >) . (a< n >,a‘">)| < e | (a<"> ,a<"> ) | (24) 

Once the induced velocity a 1 and a 2 are obtained, the 
remaining output quantities are computed by direct analytical means as 


indicated in the the following formulas: 

1 

A = 7 if A (r) dr (25) 

0 

V(r) = )[i(l a ai )] 2 + [r (1 - a 2 )] 2 ) 1/2 (26) 

M(r) = M t V ( r) (27) 

A (1 + a x ) 

<j)(r) = Arctan _ ■ ■ ^ (28) 

i 1 

T = yN / V 2 C(r) (C cos 4 > - C . sin <|>) dr (29) 

c 0 1 a 
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In 

l 


p = 

/ rCV 2 (C^ sin <f> + C d cos <(») dr 

(30) 



_0 


C T ■ 

tt 2 T / 4 


(31) 

C P = 

tt 3 P/4 


(32) 

n = 

It/p 


(33) 

y r > 

= - 

■ [r(l - a 2 ) ] 

(34) 

V z (r) 

= . 

■fA(l + a j) ] 

(35) 


Propeller Angle of Attack 

To simulate the effects of a propeller angle of attack ctp, the 
resulting flow field is rotated. This introduces an angular dependence 
of the inplane angle ip on the output quantities. Therefore, the 
thrust and inplane components of velocity become 

V^(r,i|>) = -[(r + A sin a p cos i|>) (1 - a 2 )] (36) 
V z ( r, ip) = -[A cos cx p (1 + a j) ] (37) 

Then, by using V^(r,^) and V z (r,i(»), the remaining angular 
dependent output is expressed as 


4>( r, ip) = Arctan 

(38) 

%(r,ip) 


a (r, ip) = e o + ^ - <f>( r , ip) 

(39) 

M ( r , ip) = (\ 2 + \ 2 ) 1/2 

(40) 


REFERENCES 

1. Glauert, H.: Airplane Propellers. Aerodynamic Theory, div. L, 

Volume IV, W. F. Durand, ed., Dover Publ., Inc., 1963, pp. 
169-360. 

2. Goldstein, Sydney: On the Vortex Theory of Screw Propellers. 

Proc. R. Soc. London, ser. A, vol . CXXIII, Apr. 1929, pp. 
440-465. 
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10.6 PROPELLER LOADING MODULE 


01 1 ie J. Rose 
PRC Kentron, Inc. 


INTRODUCTION 

Noise prediction for rotating blades requires a knowledge of the 
forces exerted by the fluid on the blades. These quantities are needed 
at specified points on the blade and for time intervals identified by 
basic periodicity assumptions so that a value at any space/time point 
is available by interpolation. Thus the purpose of the Propeller 
Loading Module is to produce an array of loadings at specified blade 
surface points and for specified times. 


a 

b 

C 

C 


SYMBOLS 

innermost radial station, re R 

outermost radial station, re R 

blade chord, re R 

section drag coefficient 

local friction coefficient on blade surface 

section lift coefficient 

local pressure coefficient 

speed of sound, m/s (ft/s) 

section drag force, re pR 3 ^ 2 

section lift force, re pR 3 fi 2 

net torque force acting on blade, re pR 4 n 2 , Q/Tq 

section torque force, re pR 3 n 2 

net thrust force acting on blade, re pR 4 £2 2 , T 
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F 


t 


M 


P 

Q 

R 


r 

T 

V 


V. 


* 

x,y,z 


a 

Si, 5 2 

P 

T 

$ 

* 

a 

Subscripts: 

t 


section thrust force, re pR 3 Q 2 
Mach number 
pressure, re pR 2 ^ 2 
torque, re pR 5 £2 2 
propeller disk radius 

distance of blade element from rotor axis, re R 
thrust, re pR 4 ft 2 
velocity, re Rf2 

induced axial velocity component, re Rft 

induced tangential velocity component, re Rft 

nonrotating coordinate system moving with propeller 
hub 

blade section angle of attack, rad 

blade surface coordinates 

fluid mass density, kg/m 3 (slug/ft 3 ) 

skin friction tangential stress, re pR 2 fi 2 

inflow angle, rad 

azimuth angle of blade, rad 

angular speed of rotor, rad/s 

value at blade tip 

INPUT 


The quantities required by this module can be conveniently grouped 
into three basic classes: blade geometry, flow and operating 

parameters, and aerodynamic characteristics. Blade geometry is 
furnished in the form of tables by the Blade Shape Module, the flow 
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field by the Propeller Performance Module, and aerodynamic 
characteristics by the Blade Section Aerodynamics and Blade Section 
Boundary-Layer Modules. Operating parameters are provided by the user. 


5 1 

n « 2 

c( C 1 ) 


Blade Geometry 
number of spanwise stations 
array of spanwise stations, re R 
number of chordwise stations 
array of chordwise stations, re 2 tt 
blade chord, re R 


<l» 

MUi.ilO 
“Ui.'fO 
<p(tl ,<p) 


Flow and Operating Parameters 

Do 

tip rotational Mach number, — 

c 

oo 

number of blades 

number of rotation angles 

array of rotation angles 

local effective Mach number of blade element 

local angle of attack, rad 

local inflow angle, rad 


C d (?i,«,M) 

C f Ui,5 2 ,a # M) 

C A Ui»a,M) 

C p Ui.5 2 ,a.M) 


Aerodynamic Characteristics 
section drag coefficient 
local skin friction coefficient 
section lift coefficient 
local pressure coefficient 
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OUTPUT 


The output of this module is a table of pressure and friction 
loading together with lifting-line forces. 

number of spanwise stations 

Ci array of spanwise stations, re R 

number of chordwise stations 

5 2 array of chordwise stations, re 2tt 

p(Ci>? 2 >^) pressure loading, re pR 2 ^ 2 

t(£ 1 ,Z 2 >'P) skin friction loading, re pR 2 ft 2 

N. number of rotation angles 

if) array of rotation angles, rad 

F t (C 1 ,'P) section thrust force, re pR 3 ^ 2 

Fq (C j. , ^) section torque force, re pR 3 n 2 

F T (ij)) net thrust force, re pR 4 £2 2 

Fq(<f)) net torque force, re pR 4 ft 2 

r T ( ip ) effective blade radius where net thrust force acts, 

re R 

r n (ifi) effective blade radius where net torque force acts, 

g re R 


METHOD 

The Propeller Loading Module computes aerodynamic loads on a rigid 
propeller blade, which arise as the blade operates in a given flow 
environment. Local Mach number and angle of attack are provided in 
terms of an aircraft fixed cylindrical coordinate system as shown in 
figure 1. Radial Mach number and other three-dimensional effects are 
assumed to have a negligible effect on the loading, and this allows the 
algorithm to be based on blade element theory together with 
two-dimensional aerodynamic characteristics. 
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With figures 1 to 3 as guides, equations for the resultant blade 
loads are immediately written. Consider a typical blade section of 
width dr and chord C(r) located a distance r from the rotation 
axis. The local velocity seen by this element is given by 

V(r,i|>) = M(r,i|»)/M t (1) 

From figure 3, it is evident that the thrust and torque components of 
the section loading can be written as 

F t = F i cos * " F d sin * ( 2 ) 

F q = F t sin * + F d cos * ( 3 ) 

The section lift and drag are written in terms of coefficients as 

f a = jV 2 C £ C(r) (4) 

F d = |v 2 C d C(r) (5) 

The section thrust and torque forces thus are written as 

= 2 C ( 5 1 ) {C^[^ 1 ,oi(c 1 ,ip),M(5 1 ,vp) ] cos <(> 

( 6 ) 

- C d [g 1 ,aU 1 ,i|i),MU 1 ,i|i) ] sin <j>} 

F (Ci» 4») =|v 2 0(5!) {C 4 [5 1 ,a,(Ei,i|>),M(5 1 ,i|i) ] sin 4, 

(7) 

+ C d [? 1 ,a(c 11 i|»),M(5i > i|>) ] cos 4,} 

In a similar way, equations for the surface pressure and tangential 


stress are derived as 



p( 5l > ?2> ’(’) = 

7 V 2 C p (5 1 ,5 2 ,a,M) 

(8) 

t(Ci,52»^) = 

2 ^ 2 Cf ( > ?2> a »M) 

(9) 


Net thrust and torque forces are obtained by integrating with respect 
to the radial variable. The results are given as follows: 
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( 10 ) 


b 

F t = / F t (5 lf i|i) d£ x 

a 

b 

F = / F (d.^) d^i (11) 

w a M 

Effective radii based on the net thrust and torque forces are 
calculated by the following equations: 


r T U) = -jr- 

i ' T 

V 

f 

a 

Si 

F t (?i,^) 

d Si 

(12) 

r o U) = k 

b 

/ 

a 

Si 

F q ( S i » 4>) 

dSi 

(13) 


Q a 
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Figure 2.- Blade element geometry. 





11 . PROPELLER NOISE 



11.1 SUBSONIC PROPELLER NOISE MODULE 


Sharon L. Padula 
Langley Research Center 


INTRODUCTION 

The Subsonic Propeller Noise Module computes the periodic acoustic 
pressure signature and spectrum of a propeller with subsonic tip 
speed. The computation is based on a solution of the Ffowcs 
Williams-Hawkings (FW-H) equation without the quadrupole source term. 
(See refs. 1 to 3.) The blade surface pressure and geometry are 
assumed known. The surface pressure can also be periodic in time. The 
observer is always assumed to be moving with the aircraft; only in such 
a frame is the acoustic pressure periodic. The main acoustic 
formulation uses the source distribution on the actual blade surface. 
Several approximations to the full solution are also included as 
options. 


A 


c 

d 


FFT 

f 


9 


H 


SYMBOLS 

area of blade section, m 2 (ft 2 ) 

ambient speed of sound, m/s (ft/s) 

radial distance of observer from xg-axis, m (ft) 

radial distance of source from xg-axis, m (ft) 
symbol substituted for integrands in equation (41) 
Fast Fourier Transform 
function defining blade surface 
frequency, Hz 

symbol substituted for integrands in equation (45a) 

equation g = 0 describes relationship between source 
and observer time 

azimuthal force on single blade, N (lbf) 
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k 

->■ 

L 

-► 

L* 

■¥ 

£ 


M 

m 


N 

N 

N 


b 

s 

t 


n 


P 

Ap 


PL 

Po 

Pref 

P* 


P titis 

PT 


harmonic number 

net force of blade acting on fluid, N (lbf) 

net force per unit radial distance, N/m 
( 1 bf/m) 

local force per unit area of blade acting on 
fluid, Pa (lbf/ft 2 ) 

source Mach number 

exponent; 2 m is number of time points used 
by FFT routine 

number of blades 

maximum harmonic number in calculations 

number of times 

blade surface normal vector 

surface pressure p - p , Pa (lbf/ft 2 ) 

a 0 

difference between lower and upper surface 
pressure, Pa (lbf/ft 2 ) 

absolute surface pressure. Pa (lbf/ft 2 ) 

acoustic pressure produced by loading, 

Pa (lbf/ft 2 ) 

ambient pressure. Pa (lbf/ft 2 ) 

reference pressure for definition of dB, 

20 yPa 

acoustic pressure, Pa (lbf/ft 2 ) 

total acoustic pressure produced by propeller, 

re p‘ , Pa (lbf/ft 2 ) 

K rms’ v 

root-mean-square acoustic pressure 

acoustic pressure produced by thickness. 

Pa (lbf/ft 2 ) 


r 


effective radial distance to be used in 
point source approximation, m (ft) 


distance from source 


point 


at emission time to 


observer, (r = x - y (t*) ), m (ft) 
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J + Oj '< + X + 


s 

SPL, 

k 

T 


t 


t 


v 


X 


0 

0 ' 



T 


* 


T 

4 > 

v 

<J> 


surface area 

acoustic spectrum, dB 

net thrust from single blade, N (lbf) 

time at which noise signal is received by 
observer, s 

blade surface unit tangent vector 

forward velocity of aircraft, m/s (ft/s) 

source velocity, m/s (ft/s) 

observer position in ai rcraft-f i xed frame, 
m (ft) 

observer position in ground-fixed frame, m (ft) 
source position in ground-fixed frame, m (ft) 
Dirac delta function 

source position in blade-fixed frame, m (ft) 

angle between radiation vector and surface normal 
vector, rad 

angle between radiation vector and surface tangent 
vector, rad 

root pitch change angle, rad 

elliptic surface coordinates 

ambient density, kg/m 3 (slug/ft 3 ) 

local tangential stress, N/m 2 (lbf/ft 2 ) 

time at which noise signal is emitted at source 
position, s 

solution to retarded time equation, g = 0, for 
given observer time, s 

parameter defined by equations 27 

solution to retarded time equation as solved by 
Newton's method, n ( t - t), rad 

blade volume, m 3 (ft 3 ) 

angle between x 2 - and r^-axes, rad 

angular velocity of blade, rad/s 
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Subscripts: 


L 

m 

n 

o 

r 

T 

Notation: 


★ 


component along ith coordinate axis 

(£ q -r q implies summation convention); also 
inner radius 

1 oadi ng 

property of mean blade surface 
component in direction of surface normal 
ambient condition; also outer radius 
component in direction of radiation vector 
thickness 


mean 

unit vector 
time derivative 

complex conjugate; specific solution to retarded 
time equation 


INPUT 

The computation of subsonic propeller noise requires data of the 
blade shape characteristics, aerodynamic characteristics, flow and 
operating parameters, and observer geometry. The noise signature is 
calculated by using the following methods: (1) full blade formulation, 

(2) mean surface approximation, (3) compact chord approximation, or 
(4) compact source approximation. Option numbers in parentheses 
identify the inputs which correspond to each method. 


User Parameters (1, 2, 3, 4) 

N, number of blades 

b 

N highest harmonic number desired 

R blade length measured from axis to tip, m (ft) 

V aircraft forward velocity, m/s (ft/s) 
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x + 







5i 

^1 (?1 »? 2 ) 
03(^1 >^2) 
$1 ,$ 2 

AUi) 

Pi ,£(5 i) 

P3,tUl) 


5l 


?2 

P(Ci >?2 »’!’) 
cf(5i »? 2 >^) 


P 

c 


0 


root pitch change angle, rad 
initial azimuth angle of first blade, rad 
angular velocity of blade, rad/s 
blade volume, re R 3 

Blade Shape Table (1, 2) 
spanwise location, re R 
chordwise location, rad 
blade surface abscissa, re R 
blade surface ordinate, re R 
slopes of basis functions at nodes 

Spanwise Function Table (3) 
blade section area, re R 2 

abscissa of leading edge, re R 
ordinate of leading edge, re R 

Aerodynamic Characteristics Table (1, 2) 
spanwise location, re R 
chordwise location, rad 
azimuthal location, rad 
pressure loading, re p 0 (Rft) 2 
skin friction loading, re p 0 (Rft) 2 

Observer Table (1, 2, 3, 4) 
initial observer position, m (ft) 

Atmospheric Properties Table (1, 2, 3, 4) 
ambient fluid density, kg/m 3 (slugs/ft 3 ) 
local sound speed, m/s (ft/s) 
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F z ( ? i , 4») 

T( t|i) 

H(t|») 

5 t U) 

EqU) 


Performance Characteristics Table 
axial section force, re p Q R 3 fl 2 (3) 
azimuthal section force, re p Q R 3 fi 2 (3) 
net thrust for a single blade, re p Q R 4 fi 2 (4) 

net azimuthal force for a single blade, 

re p R 4 ft 2 (4) 
o 

effective radius based on net thrust, re R (4) 
effective radius based on azimuthal force, re R (4) 


OUTPUT 


This module provides the complex pressure spectrum and pressure 
time history at a fixed observer location measured . from the center of 
the hub. Complex pressures are given with the e“ int time harmonic 
convention and all spectra are understood to be two sided with 
p (x,-k) = p*(x,k). With this convention, the mean-squared pressure is 


2p p* 

P P 


for 


each 


harmonic (where 


2p'p'* = P 
PP 


rms 


). 


k 

f k 

c k (>0 

SPL k (x) 


Complex Pressure Spectrum 
harmonic number 
frequency, Hz 

complex acoustic pressure, re p 0 (Rn) 2 
sound pressure level associated with C, (x) 


Pressure Time History 
t time, s 

p p (x,t) acoustic pressure, re p rm$ 
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METHOD 


Acoustic Formulation 

The governing acoustic equation is the following inhomogeneous 
wave equation: 


L- : iV 
c 2 Vt 2 


’V = !t tVn| 7f | S(f> ] 


" [*i| 7f | s ( f )] 0) 

This is the FW-H equation without^the quadrupole source term. The 
blade surface is described by f(x,t) = 0. The acoustic pressure is 
denoted by p', v n is the local normal velocity of the blade, and 
£j is the local force per unit area of the blade acting on the 
fluid. The density and the speed of sound in the undisturbed medium 
are denoted as p 0 and c, respectively. The symbol 6(f) stands 
for the Dirac delta function. The two terms on the right of 
equation (1) are known as the thickness and loading noise sources, 
respectively. Note that equation (1) is written in a Cartesian frame 
fixed to the undisturbed medium which is called the ground-fixed 
frame . Two other frames of reference are more suitable for describing 
the solution and computation of the noise. 

The FW-H equation is valid in the entire three-dimensional 
unbounded space. The Green's function of the wave equation for the 
unbounded space is 6(g)/4irr where g = t - t + r/c and r = | x - y | . 

Here i and t are the source and the observer times and y and x 

are the source and observer positions, respectively. With this Green's 
function, the formal solution of equation (1) is 

4irp'(x,t) = 1 p 0 v n | vf 1 6 ( f ) «(g) d y dT 

'it S? *i| vf l 6(f) 6(g) dT (2) 

Note that the integrals in equation (2) are four-dimensional as written 
bel ow: 

/ ... dy dr = f- I" f°° f°° ... dy dr (3) 

.00 — OO .00 .00 


Before integrating equation (2), the following identity, 
established by differentiation, is used 



r$(g)l _ 1 9 

r\j<5(g) 

|_ r J c at 

r 


r i 6(g) 


+ 


(4) 


where r-j = (x-j - y-j)/r is the unit vector in radiation direction. 

Since only the function 6(g)/r is dependent on x, the second term in 
equation (2) may be rewritten as 


3x 


- / 7 ^ vf 6(f) 6(g) dy d T = | / 7 V v f «( f ) «(g) d y d T 


+ / - Z vf 6(f) 6(g) dy d T (5) 
r 2 r 

A 

where Z r = t-jr-j. Using equation (5) in equation (2) gives 

4 tt P 1 ( x ,t ) It / 7 (p 0 cv n + v ) | Vf | 6(f) 5(g) dT 

+ / ~ Z r | vf | 6(f) 6(g) dy d T (6) 


Now a Cartesian frame, denoted as the r[-frame, fixed to the blade 
is introduced. This frame is referred to as the blade-fixed frame . It 
is assumed that the blade is rigid so that its surface description in 
the ri-frame is time independent. To integrate equation (6), the 
following steps are employed: 

(1) Use the transformation of variable y n 

(2) Use the transformation t -*■ g to integrate with respect to 

the variable g 

(3) Integrate the remaining delta functions 

Step 1 : Since the Jacobian of the transformati on is 1, it follows 

that the volume elements in the two spaces are equal. Also note that 
in the ri-frame the equation of the blade surface depends on rf only. 
This is obtained as follows: 

f [y(n,T),x] = f(n) = 0 (7) 

That is, the equation of the blade surface in the n-frame is obtained 
by substituting for y in terms of n and t. Note that since the 

blade-fixed frame is in motion, the transformation y -*• rj is time 

dependent and y = y(n,x). To reduce confusion, the tilde on f ( n) is 
dropped in the following results. The function g is defined as 
follows: 


g = x - t + 


X 


y(n,T) 


/c 


( 8 ) 
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Step 2: Now fix n and use the transformation x ■* g. The 

Jacobian of the transformation is 


J 




( 9 ) 


where M r = v-jr-j/c. Here v is the velocity of the point with 
position vector n in the blade-fixed frame relative to the 
ground-fixed frame. In step 3, this velocity will be defined more 
carefully. The result of steps 1 and 2 is 


4xrp ' (x 


t> .ii./ 
c at J 


P„CV + 

o n r 
r 1 1 - M. 


| vf | 6(f) dn 


/ 


r 2 1 - 


r|Jx 


Vf 6(f) dn 


( 10 ) 


Here the integrals are evaluated in the blade-fixed frame over the 
entire space. The source time x* is the emission time obtained by 
finding the root of the equation 




g = t 


t + 


y(n,t*) 


/c = 0 


( 11 ) 


For subsonic propellers considered herein, this equation has only one 
root. 


Step 3: In this step, the integrations over the delta functions 

in equati on (10) are performed by using the relation 


df dS 
i vf j 


( 12 ) 


where at this stage dS is the element of surface area of the surface 
f = Constant. This relation is derived in reference 2. Substituting 
equation (12) into equation (10) and integrating with respect to f 
result in the following equation: 


4irp'(x,t) 


I i_ 

c at 



p 

cv + i 

i 


— 

£ 

f 


) n 

r 

i . . 

dS + / 
x* f =0 


r 

J 

=0 

r 

1 - M 

r 

r 

r 2 

11 - M 

r 


dS (13) 


Here, M r = v-jr-j/c, where v is the local velocity of points on 
the blade surface itself. Since M r < 1 for all points on a 
subsonic blade, the absolute value signs about 1 - M r are dropped in 
the following formulas. 
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Equation (13) has been used for high-speed propeller noise 
calculations by Farassat. (See ref. 4.) Since for subsonic propellers 
M r < 1, the time derivative can be taken inside the integral. The 
following terms are functions of time in the first integral of equation 
(13): 


r( 1 - M p ) 


r(l - M 


The chain rule of differentiation is then used to differentiate the 
integrand of the first integral in equation (13): 

9 _ = _ 3 _ . 

at at at* u 

From equation (11), by differentiating with respect to t, the 
following is obtained: 


1 - M • 
r 


Therefore, from equation (14), 


at l - M ax* 

r 


The following relations are thus established: 


9r i _ r i v r - v i 
3x* r 


^r _ l_ r + 2 21 

3x* cr f r i 3x* v r " v ] 


By using equations (16) through (19), the time derivative in 
equation (13) is carried out explicitly to give: 
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4 Trp. 1 (x,t) = / 


f=0 


0 . r . 
*1 1 


r( 1 - M p ) 2 


dS 


* / 

f=0 


V - £ i M i 
r 2 (l - M p ) 2 ^ 


dS 


+ I 

c 


/ 

f=0 


£ (rM. r. + cM -c M 2 ) 
r v ii r _ 

r 2 (1 - M ) 3 
r 


dS (20a) 


and 


47TP-J. ( X , t ) = 


/ 

f=0 


p o v n( rM i^i + cM r " cN|2 ) 

r 2 (l - M ) 3 
r' 


dS 


P * ( X, t ) = p^(x,t) + P-J-(x.t) 


(20b) 


(20c) 


where p^ and p| denote the loading and thickness noise, 
respectively. Equations (20) are used when the sources on the 
full-blade surface are specified. They are related to an equation 
derived by Woan and Gregorek. (See refs. 4 and 5.) 

The symbols and are defined as follows: 


v . 

M. = — 
1 c 


(21) 


v . 

• l 

M. = — 
1 c 


( 22 ) 


Note that v, and v i are the local velocity and acceleration of the 
blade with respect to the ground-fixed frame. It is assumed here that 
the propeller blades are rigid. It is important to note that £-j is 
the rate of change of the vector £-j as seen from the ground-fixed 
frame. This means that even if is steady in the blade-fixed 

frame, £j is not zero. 

The force intensity vector £-j can be written as 

£-j = P n i + crt. ( 23 ) 

where p is tie surface pressure defined as p a - p 0 . Here p a 
is the absolute surface pressure and p 0 is the ambient pressure in 
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the undisturbed medium. The local tangential stress is denoted by a 
and the unit tangent vector pointing toward the section leading edge is 
denoted by £,. From equation (23), the rate of change of is 

found to be 



£. = pn i + pn i + at. + at i 

(24) 

shown 

that 


• 

a 

n . = 

l 

(fi x n) i 

(25a) 

• 

l i = 

(fi x t)^ 

(25b) 


From equations (25) and (24), the following results are obtained: 


£ r^ = p cos 0 + a cos 0'+ pT x + aT 2 (26a) 

r.j = p cos 0 + a cos 0 1 (26b) 

£.M. = pM + oM. (26c) 

li r n t 


Note that 0 is the angle between n and r, o' is the angle 
between £. and r. , and T 1 and T 2 are defined as follows: 


i = r • ( 

ft x n) 

(27a) 

2 = r • ( 

< -4— > 
X 
+CZ 

(27b) 


Also M-t is defined as v-jt-j/c. Equations (26) are substituted 
into equations (20) and the result is solved numerically. 


A useful approximation to equations (20) is written next. If the 
propeller blades are assumed to be thin, the acoustic sources can be 
moved to the mean surface of the blades. The mean surface is defined 
to be that surface formed by generating the section mean curves along 
the leading-edge curve. Let the surface be denoted by the equation 

f* (x,t) = 0. The FW-H equation becomes 
m v ’ ' 


1_ 9 2 p 1 

C 2 3t 2 


V 2 P ' 



ml 

vf 


m 


I 5 < f m >] 


(28) 
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where Ap = Piower “ Pupper anc * n im is the unit normal to 
f = 0 pointing toward tne upper side of the blade. Here v n is the 
1 Seal normal velocity of the blade on one side of the mean surface due 
to thicknes_s distribution alone (i.e., the camber effect is ignored). 
Note that v n is calculated by evaluating v n at each upper 
surface point and v n at the correspondi ng point on the lower 
surface and then averaging. Furthermore, skin friction is neglected. 


Following a procedure similar to that leading to equations (20), 
an expression for p‘(x,t) is 


4irp^(x,t) = 


i J 


f =0 
m 


Ap cos 0 + ApT, 
r(l - M r ) 2 


dS 


f 

Ap ( COS 0 - M n ) 

J 

f =0 

r 2 

(1 - M p ) 2 _ 

m 


I 

r 

Ap COS 0 ( 


f =0 
m 


dS 




r 2 (1 - M r ) 3 


dS (29a) 


4irpj(x,t) 


/ 


f =0 
m 


2 p 0 v n (rM.r i + cM p - cM 2 ) 


dS 


(29b) 


r 2 (l - M ) 3 
r ' 

p'(x,t) = P^(x,t) + p-j-(x.t) (29c) 

In equation (29a), 0 is the angle between the radiation vector r 


and n. . 

i m 


Two compact-source approximations for the subsonic propeller noise 
calculation are also used. In the first approximation, each blade is 
replaced by a single point source involving both the thickness and 
loading effects. The loading noise formula of Lowson and the thickness 
noise formula of Succi are used. (See refs. 6 and 7.) Let L, be 
the net force by each blade on the fluid and v be the net volume of 
each blade. Then the loading and thickness noise are 
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4-rrp^ (x ,t ) = 


icr(l - M p ) 2 


? i L V irh r )[Vr (1 - m2 h - 7 M i L i 


» F L (x,t,R E ) 


and 


p v 

4irp|(x,t) = ' 


(1 - M r )( M p - 3 \ M-M) + 3(p- M p ) 2 


r(l - M ) 5 L 

+ 3(M p - j M 2 ) 2 + £ (M p - j M 2 )(l + 4M f + M 2 ) 
= F T (x,t,R E ) 


T=T* 


T= T* 


(30) 


(31) 


Here the sources are assumed to be located at some effective radial 
distance (usually about 0.8 radius of the propeller). The effect 

of each blade is computed separately and summed to get the total noise. 

A second compact-source formulation is obtained by assuming that 
each blade can be approximated by a line source (chordwise compactness 
assumption). In equations (30) and (31), let LMCj) and A(? : ) be 
substituted for L. and y, respectively, and represent the force per 
unit distance and ilade volume per unit distance. Here is the 

radial distance in the blade-fixed frame. The noise due to loading and 
thickness of the propeller is then given by 


4irp'(x,t) = / 0 [F. 1 (x,t,Cj) + Fi(x,t,?i)] dc x 
R i 

where the integrands F^ and Fj are defined by 


(32a) 


FMx.t.gj) = 


cr(l - M ) 2 
r 


l r Vi + 


r . L 

1 1 [M p + j (1 - M 2 )] 


(1 - M p ) 


? M i L i> 


(32b) 
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and 


Fy (x,t, 5l ) 


p 0 A(5i) 
r(l - M p ) 5 



M )(M 
r /\ r 


3 - M-M) + 3(— M ) 2 
r r r' 


+ 3 (M p - -pM 2 ) 2 + -p (M p - ± M 2 ) ( 1 + 4M p + M*) 


(32c) 


The symbols R-j and R 0 in equation (32a) represent the inner and 
the outer radii of the blade. The effect of several blades is taken 
into account by a method requiring a single blade calculation as 
discussed later,. 


Frames of Reference 

In this section, the coordinate systems used in the acoustic 
calculations will be discussed. Basically there are two reference 
frames: the ground-fixed x-frame and the blade-fixed n-frame. The 

x-frame which remains fixed to the undisturbed medium is set up as 
follows: 


The origin of this frame is at the propeller center at time 
t = 0 

The propeller disk is in the x 1 x 2 -plane and the x 3 -axis 

coincides with the propeller axis, with positive direction in 
the flight direction 

The xj-axis is assumed to be orthogonal to the x 3 -axis, upward 
and parallel to the plane of symmetry of the aircraft 

The x 2 -axis is defined in such a way that the x-frame is 
right-handed (see fig. 1) 

The observer position is always specified in this frame as is explained 
shortly. 

The blade-fixed frame is set up as follows: 


The origin coincides with the origin of the x-frame at t = 0 

The n 3 -axis coincides with the x 3 -axis 

The n 2 -axis lies on the pitch change axis of the blade 

The mn 2 -plane is normal to the propeller shaft 

The n x -axis is defined to make the n-frame right-handed (see 
fig. 2) 
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Note that the blade mean surface does not lie in the n^-plane but 
makes an angle equal to the geometric pitch angle at each radial 
position. Note also that the in- frame rotates with angular velocity 0. 

In the acoustic formulas presented in the previous section, many 
operations involving dot products or cross products of two vectors are 
required. When the products involve blade normals or surface 
pressures, it is convenient to transform the components of these 
vectors to the ri-frame and then perform the operation. However, note 
that the source position y always refers to a source point in the 
x-frame with a known position ri in the blade-fixed frame. The 
quantities r = x - y and r = r/r are calculated in the x-frame, and 
then the vector components are calculated in the ri-frame. Note that 
for any variable which has a velocity term, such as M r = M-j r-j 
and v n = v*n, the velocity is with respect to the x-frame. In other 
words 


v i 


*1 


(33) 


The components of this vector must then be calculated in the ri-frame to 
get M r and v n . 

Another frame of reference which is used conceptually is the 
ai rcraft-fi xed nonrotating )t-frame. This frame coincides with the 
ground-fixed x-frame at t = 0. The acoustic pressure p'(£,t) for 
an observer in the aircraft-fixed frame with the position vector )T is 
found as follows: Since both )t and t are known, the relative 

position of the ground-fixed frame and the ai rcraft-fi xed frame is 
known. From this, calculate the observer position x in the 
ground-fixed frame. By using this x and t, calculate p‘(x,t) 
which is the same as the acoustic pressure at position t in the 
aircraft-fixed frame and at time t. Thus, the )t-frame itself is never 
used in acoustic calculations, even though the pressure signal which 
results is p 1 ()t,t) . 


The transformation from one frame to the other at time t, in 
matrix notation is 


~*r 


cos 4) -sin ip 0 


’ll 

x 2 

= 

sin ip cos <Jj 0 


n 2 

^3 


0 0 1 


^3 

— 


— 





(34a) 
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A final coordinate transformation is necessary to account for 
changes in the root pitch setting. This rotation through a small 
angle e 0 is applied to the blade surface coordinates, which are 
calculated by the Blade Shape Module: 



Whenever the n-frame or pitch axis coordinate system is mentioned 
herein, it is assumed that the above corrective transformation has 
occurred. Often the Blade Shape Module and Subsonic Propeller Noise 
Module will be executed without change to the collective pitch 
setting. Thus, e 0 = 0 and the last transformati on is unnecessary. 


Computational Strategy 

A pressure time history of a multibladed propeller is calculated 
by summing the acoustic pressures of each blade as a function of time. 
The signals are periodic because the observer position is a fixed 
distance from the hub at all times. Thus, the signal produced by 1 
revolution of a single blade can be used to determine the time history 
and spectrum of the propeller. The general computational strategy is 
outlined below. Details concerning the different methods for 
calculating p'(t,t) are contained in later sections. 

The steps involved in calculating the acoustic time signature for 
a single blade are as follows. Assume that the initial pitch change 
axis of the blade is aligned with the X 2 -axis. That is, assume 
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t = ip = 0 and the x-frame, X-frame, and n-frame coincide. Divide the 
time required for 1 full revolution of the blade into even time 
increments such that 


= (j ~ *) 27r 

i 


(i = 1, 2, ..., N t ) 


Likewise, represent the blade by a set of source points n(j), where 
the number and location of the source points depend on which method for 
calculating noise is selected. Then, for each time t-j and^each 
source point n(j), determine the source y and observer x 
locations in the ground-fixed frame. With these, solve the retarded 
time equation and calculate distances and velocities at emission 
time. Next, calculate vector components in the rotating n-frame and 
calculate tiie influence of this source point on the total noise. 
Finally, accumulate the information for each time and source point into 
a complete time history. 

The steps involved in producing a time history and spectrum for 
the propeller are as follows: 

Tabulate the time history for one blade so that intermediate 
values can be obtained by interpolation. 

Decide how many evenly spaced time points are required by the 
Fast Fourier Transform (FFT) procedure. If N s is the 
highest harmonic number of interest and N^ the number of 
time points, find an integer m such that 2™ is at least 
as large as N^ _> 2N S . 

Calculate the total acoustic pressure time history by adding 
the contribution from each blade. 


Finally, perform an FFT to get the pressure spectrum. 


The above process can be stated mathematically as follows. The 
time required for 1 revolution of the blade is divided equally into 
2 m points such that 


= (k - 1 ) 2 TT 

k (2 m )fi 


(k = 1, 2, 



m = INT (log 2 N t ) + 1 


(36) 


where the INT function turns a real number into the largest integer not 
exceeding the real number. 
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The acoustic signal generated by a propeller is just the sum of 
the pressures produced by each blade. The pressure produced at time t 
by blade b, which had an initial angle tfo = + 2ir(b - 1 ) / 

relative to the X 2 axis, can be found by interpolating into the single 
blade pressure table for p 1 (jt.t+ipb/fi) • Thus, the total acoustic 
signal is 

N h 

> k -»■ m 

P T (X,t k ) = y p'(X,t k + (k = 1 , 2, ..., 2) (37) 

b=l 

-V 

where X is the observer position in the aircraft fixed frame, is 

the number of blades, and ^ is the initial position of each 
blade. The spectrum associated with this time history is calculated by 
using a standard FFT routine. Appendix A contains further details. 


Retarded Time Equation 

This section describes a method for solving the retarded time 
equation, g = t - t + r/c = 0. As noted previously, the root t = t* 
of this equation must be found for each source point used to represent 
the blade and for each observer time t. Since the roots of the 
equation are calculated repeatedly it is imperative to solve the 
equation in an efficient way. Details of the solution technique are 
reserved for appendix B. 

The problem can be stated in terms of arbitrary receiver location 
x and source location y defined in the ground-fixed coordinate system 
at time t. As described in reference 4, the equation is written 

r 2 :: x 2 + d 2 + d 2 - 2d d cos (fit + \i> - ip ) 

3 x y x y v y y x 

+ V f 2 t 2 - 2x 3 Vp t 

= : c 2 ( t - t) 2 (38) 

where d x and dy are the radial distances between the x 3 -axis 
and the observer and source points, and ii> x and jjy are the 
observer and source angles of rotation at t = 0. For instance, if the 
source point is located on the pitch change axis then dy = n 2 and 
= 0. (See fig. 2.) After some manipulation, 

r 2 = ( x 3 - V p t) 2 + d 2 + d^ - 2d x d y cos [n( T - t) + ^ - (^ - nt) ] 

+ V 2 ( x - t) 2 - 2Vp (x 3 - Vp t ) ( t - t) 

= c 2 (t - t) 2 (39) 
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Letting <}> = fi{ x - t), equation (39) is written as 


A<f > 2 + B 4 > + C + cos ( <j> + D) - 0 


( 40 ) 


where 


A = (c 2 - V 2 )/2fi 2 d x d y (A > 0) 

B = V p (x 3 - V p t)/^d x d y 
C = ~ [ ( x 3 - V p t) 2 + d 2 + d 2 ]/2d x d y 

D = lb - lb + fit 


Based on the transformation given in equation (34b), the term 
(x 3 - Vpt) can be replaced by the constant X 3 for all times. 
Consequently, the coefficients A, B, and C in equation (40) change only 
if source position changes. Thus, the solution to equation (40) for 
one source point and source time will be a good initial guess for the 
same source point at a slightly later time. 


Equation (40) can be solved very efficiently by using Newton's 
method. Only one root such that (t - t) < 0 is physically possible 
in the subsonic case. With the initial guess outlined in appendix B, 
Newton's method will converge to the proper root. 


Full Blade Formulation 

This section discusses the calculation of the periodic acoustic 
pressure signal for a single rotating blade with subsonic tip speed. 
The method described requires a detailed knowledge of the blade 
geometry and blade surface pressures. The three sections which follow 
discuss simplifying approximations to the full blade formulation which 
are used in the absence of full blade information. 

The acoustic pressure at time t for an observer at point x is 
given by equations ( 20 ) as 
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4irp^(x,t) = 


/ 

f= 0 


Vi 


r(l - M ) 
r ' 


dS 


+ / 

f=0 


V - Vi 

r 2 (1 - M ) 2 


r Jt 


dS 


+ r / 


f=0 


Z (rM.r. + cM - cM 2 ) 
r v i i r _ 

r 2 (l - M ) 3 
r 


dS (41a) 


4Trp'(x,t) = / 

1 f=0 


p v ( rM . r . + cM - cM 2 ) 
n v in r 

r 2 (l - M p ) 3 


dS 


(41b) 


with 


p'(x,t) = p£(x,t) + p j ( x , t ) (41c) 

Using symbols Fl \ , Fl 2 , and F[_ 3 to stand for the integrands 
which are evaluated at retarded time x*, allows equations (41) to be 
written in short form as 


4ir Pi = 7 I F|_i( h) dS + / Fi_ 2 ( n) dS + ~ 

L c f=0 f=0 c 


/ F L3 (n) dS 
f=0 


(42a) 


4 TTp-1 = / F T (n) dS (42b) 

1 f=0 1 

Using the Blade Shape Module, the surface of a propeller blade is 
described by a set of points n(j). The three components of each point 
n(j) are tabulated with respect to spanwise and chordwise station 
parameters and ^ 2 . (See fig. 3.) Thus each term in equations 
(42) may be replaced by an equivalent double integral over and 

5 2 . For instance, 


/ F Ll ( n) dS = / / F Ll (n(?-i,5 2 ))J d ?2 d^ 

f=0 0 0 

where J is the Jacobian of the transformation. 


J 


3n x 3n 

3^1 3£ 2 


(43) 


(44) 
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Equation (43) is evaluated by using a compound rule. Thus, 


1 2 IT 1 

/ / F Ll J d? 2 d£ x = / Gl 1 (c 1 ) d£ x (45) 

0 0 0 

where 

2 7T 

g Li(?i) = / f Li(Si>? 2 ) J(5i. 5 2 ) d ?2 (46) 

0 ■ 

Each of the integrals over a single parameter is evaluated by using a 
standard procedure based on cubic splines. 

^ Conceptually, the process of evaluating equations (41) to get 
p*(x,t) involves calculating the values of each integrand at a set of 
blade surface points n(j) and accumulating them in a way which gives 
a sufficiently good approximation to the integral. 


Mean Surface Approximation 

In this section, approximation formulas for noise calculations 
based on the assumption of a thin airfoil are presented. The mean 
blade surface is defined to be that surface generated by the mean 
section curves. The calculation procedure is derived from 
equations (29) which are 


4Trp^(x,t) = - / 


f =0 
m 


Ap cos 9 + Ap T 1 
r ( 1 - M p ) 2 


dS 


/ 


f =0 
m 


Ap ( COS 0 - M ) 

r 2 (l - M p ) 2 


dS 


/ 


f =0 
m 


Ap cos e (rM^r.. + cM^ - cM 2 ) 


r 2 (1 - M ) 

v r / 


dS (47a) 


4TCPy(x ,t ) = / 


f =0 
m 


2 p 0 v n ( r Vi + c \ “ cM2 ) 


r 2 (l - M p ) : 


dS 


(47b) 


p'(x,t) = p^(x,t) + p|(x,t) 


(47c) 


Equations for the various terms of the integrands as functions of known 
quantities are determined as follows: 
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Coordinates of mean surface: 

For a given radial distance, ? ]5 f 2 varies in the range 
0 ± ^ and de fi nes a mean section line as 

'n-i ( 5i ,C 2 ) =7 [o i (C 1 ,5 2 ) + n i (5 1 ,2 1 r-5 2 ) 1 (i = 1, 2, 3) (48) 

Differential element on mean surface: 

In terms of the independent variables e >1 and _ ? 2 , the differential 
area element dS is given as 

dS = JU lf £ 2 ) d^ d ?2 = /EG - F 2 d El d^ 2 (49) 

The symbols E, F, and G are written in terms of the mean surface 
coordinates as 


® /9n 2 \ 2 

yti) + \^ij 


9 ri 3 V 


(50) 


3 Ti ^ 9^2 9^2 9^3 3^3 

F = + + 9^7 'Wl 


(51) 


G = 



9V 2 

9C 2 



9n 3 N2 


(52) 


Since n 2 = Sj and n 2 is independent of £ 2 , equations (50) and (51) 
simplify to 


/ — \ r 


‘ * m * ‘ 


(53) 


and 


9t"I 1 9n i 9 113 9n 3 

F = 9?7 W 2 + W7 W 2 


(54) 


Velocities and Mach numbers: 


The position of a point on the blade is obtained by using the 
transformation defined by equations (34) 
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cos fit -sin fit 0 


sin fit cos fit 0 


The velocity of the point is 


-sin fit -cos fit 0 


y 2 = fi cos fit -sin fit 0 


The position of a source point on the mean surface is, by analogy with 
equation (55), 


cos fit -sin fit 0 


sin fit cos fit 0 


The unit normal to the mean surface is 


and the unit vector in the radiation direction from a point on the mean 
surface is written as 


x - y 


x - y 


By analogy with equation (56), the velocity of a point on the mean 
surface is 
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-sin fit -cos fit 0 


v' = fi cos fit -sin fit 0 



The acceleration therefore becomes 


(It “ 


-cos fit sin fit 0 


-sin fit -cos fit 0 



Hence the local Mach number and its time rate of change are 


-sin fit - cos fit 0 


— cos fit - sin fit 0 
c 


n 2 + c 0 



cos fit sin fit 0 


-sin fit -cos fit 0 



The mean normal velocity and radiation Mach number are written as 


'n = 7 


[v(?!,C 2 ) * n(? 1 ,? 2 ) + v(5 lf 2ir - ? 2 ) • nU^* - ? 2 ) ] (64) 


M = M 
r 
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Miscellaneous integrand terms: 

The differential pressure Ap is defined as 

Ap - pUi,C 2 ) " P ( C i , 2 TT- 5 2 ) ( 66 ) 

The radiation angle is given by 


cos 0 = n • r 

and the symbol T x is defined to be 


(67) 




(a x n) • (x - y)/ 


x - y 


( 68 ) 


In summary, equations (47) are expressed as double integrals in the 
variables £, l and £ 2 . The equations are 


TT >•! 


4ttp' ( x , t ) = - ~ f f 


0 0 


Ap COS 0 + Ap T t 
r( 1 - M ) 2 


J(Ci,C 2 ) d ? 2 


T A 


- f f 

0 0 


’Ap(cos 0 - M )' 
r n_ 

_ r 2 (1 - M p ) 2 


J(Ci,S 2 ) dCj d ? 2 


,ir -1 


/ / 

0 0 


Ap cos 0 (rM.r. + cM - cM 2 ) 
r i i r 


r 2 (1 - M p ) 


J(?i,C 2 ) dgj dc 2 


I T 


(69) 


and 


n A 


4 irPy ( X , t ) = f f 
1 0 0 


2 p v ( rM . r . + cM + cM 2 ) 
M o n l l r 


r 2 (1 - M ) 3 
r 


J ( 5 1> ? 2 ) d? x d? 2 


—I T 


(70) 


where J = /EG - F 2 , with E, G, and F having the definitions given 
by equations (50) to (51). 


Compact Chord Approximation 

In this section, approximation formulas based on the assumption of 
a line source are presented. The basis for the calculation procedure 
is the set of equations found in equations (32): 
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4-TTp ‘ ( X , t ) = 


(71) 


/ ° [F‘ (x.t.Sj) + F‘ (x,t,Ci) 1 d^ 
R i 


where 


F^x.t.Sj) = : 


1 




r L! + [M + f (1 - M z ) ] 

cr (1 - M ) 2 1 (1 - M ) 


r i l 


(72) 


and 


F|( x ,t , c ! ) = 


P-A^j) f .. r + T r 

0 1 (1 - M )( M - 3 7 - M *M) + 3(£ M ) 2 


r(l - M r ) 5 


r' r 


r r 


+ 3(M - -M 2 ) 2 + - (M - - M 2 )(l + 4M + M 2 ) \ (73) 

v r p ' r v r r r r'J 


->■ -V 

In equation (71), the quantities M, M, M r , and r are computed as 
in equations (55) to (64) of the mean-surface approximation section. 

The leading-edge curve of the propeller is taken as the line source, so 
that all positions r are on the leading-edge curve. The quantity 
A(? x ) is the cross-sectional area of a given blade section. 

+ The key quantities to be calculated in equation (72) are l 1 and 
dL'/dt, which represent the loading per unit length and its time rate 
of change. The loading per unit length is given by 


~i ~~ 


cos fit F^(£ ls t) 

L 2 

- 

sin fit F Ujjt) 

L 3 


-F z Ui,t) 

. 


_ 


(74) 


where and F z are provided by the Propeller Loading Module. 

The time rate of change of L‘ becomes 
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*” ~ ~ 


3F 


Li 


-fi sin nt F, + cos fit — 

ip at 




aF 

d 

dt 


= 

+fi cos fit F, + sin fit —rr- 
ip at 




aF 


i 1 


z 


^ 3 | 


at 


( 75 ) 


The quantities 


3 F 8t and 


3F Z / 3t must be calculated as needed. 


Point Source Approximation 

In this section, approximation formulas based on the assumption of 
a single noise source representing each blade are presented. The basis 
for the calculation procedure is equations (30) and (31) which are 


’ { ^r-T T T 7 d - " 2 >] 


- M.L. 

r T = T ' 


(76) 


and 


4irpj(x,t) = 


Po" 


r(l - M p ) : 


(1 - M p ) (M p - 3 £ M.M) + 3(— M p ) 2 


+ 3(M - — M 2 ) 2 + — (M - — M 2 ) ( 1 + 4M + M 2 ) 
v r r r v r r M r r' 


(77) 


T=T* 


The local quantities such as r, M, and M r are evaluated at some 
effective radial distance, = R^, usually about 0.8 of the 
propeller radius (appendix C). The effective blade load is given by 


Li 


H cos fit 

L 2 

= 

H sin fit 

L 3 


-T 

___ - 





( 78 ) 


11.1-28 


where H and T are the azimuthal force and net thrust from a single 
blade and are provided by the Propeller Loading Module. The time 
rate of change of C becomes 



• • 

The quantities H and T must be calculated as needed. If it is 
assumed that the effective source point is located on the n 2 -a*is at 
the effective radius , then the local Mach number at this point 
i s gi ven by 



The first and second time derivatives are 
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The position of the source is 


yi 


cos fit -sin fit 0 


0 


0 

Yl 

= 

sin fit cos fit 0 


Re 

+ 

0 

y-i 


0 0 1 


0 


V F t 


( 83 ) 


and the unit vector in the radiation direction is given by the 
expression 



x - y(R E ) 
X - y(R E ) 


(84) 
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APPENDIX A 


ACOUSTIC SPECTRUM CALCULATION 


The acoustic spectrum is calculated from the amplitude of the 
complex coefficients of the Fourier series for pp(x,t). A standard 
Fast Fourier Transform routine based on the Cooley-Tukey algorithm is 
used for this purpose. 

Let N s be the number of harmonics desired. Find m such that 
2 m-l < 2 N s <_ 2 m and form the complex array: 

Zj = (Pp(x,t j ) ,0) (j = 1, 2, .... N t -1) (Al) 


where 


N t = 2 m t. = (2Trj/N t ) 

Now find the inverse Fourier transform 

N -1 

! t - 2 Tri kj'/N 

C. = l Z.e (k = 0, 1 N ) 

K \ j = 0 j t 

The acoustic spectrum is given by 

SPL k = 10 log (2|C k | 2 /p2 ef ) (k = 1, 2 N $ ) 


(A2) 


(A3) 


(A4) 
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APPENDIX 8 


ROOTS OF RETARDED TIME EQUATION 


The retarded time equation is written as 

A<j> 2 + B<f> + C + cos (<}> + D) = 0 (Bl) 


where only the root U/fi) <0 is of interest. Since both 
equation (Bl) and its first derivative have an oscillating term, a 
close initial guess is required to ensure the convergence of Newton's 
method. This guess is accomplished by replacing the cosine term by a 
parabola and solving the resulting quadratic equation. Details are 
given as follows. 


Step 1 - Bracket the root : 

Figure Bl illustrates that the roots of equation (Bl) are the 
intersection of a parabola and a cosine curve. As a first 
approximation to the root, find <j>j such that 

A<t>! 2 + B*! + C = 0 ( <f> x < 0 ) (B2) 


Next determine which half'cycle of the cosine curve to replace by a 
quadratic. 

Let 

f INT [Uj + D)/tt] (<h + D) >_ 0) \ 

l = < / (B3) 

INT [( <*>! + D ) / tt ] - 1 U 2 + D) < 0) J 

where Z is the integer half-cycle number and INT stands for the 
truncation integer function. 


Step 2 - Determine initial guess : 


The cosine term is replaced either by a straight line or by a 
parabola depending on whether z is an even integer or an odd integer 


cos (<)> + D) 


- <p - D + ( Z + 1/2)tt 
1 - 4[ <f> + D + ( £, + 1) n] 2 / it 2 


( z even) 
(z odd) 


(B4) 


With this substitution, equation (Bl) is solved by using the quadratic 
formula to give an initial guess <j> 2 - 
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Step 3 - Newton's method : 


The roots of equation (Bl) are determined through an iterative 
process. 



[A^ + 8<j> n + C + cos (<fr n + 0) ] 
[2A$ n + B - sin (<|> n + D) ] 


(B5) 


which terminates wnen <j> n +i - <j> n | is smaller than a user- 
specified convergence criterion or if the number of iterations exceeds 
a limit set by the user. 


Step 4 - Return answer : 

If Un+l/ fi ) >0 or if <(> n+ i substituted in equation (Bl) 
does not return a value less in magnitude than a user-specified 
tolerance, then an error condition results; otherwise, 

<f> n +i = - t). 

Note that the retarded time equation is solved for one source 
point at many receiver times. Thus, 4 > n+ i may be a close initial 
guess to the next required root. In this case. Steps 1 and 2 are 
omitted. 
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APPENDIX C 


EFFECTIVE POINT SOURCE RADIUS 


Let T be the net thrust for a single blade located at radius R_ 
and H the net inplane force located at R., (fig. Cl). A resultant 
force Fp and corresponding effective radius R F are calculated based 
on the concepts of equivalent forces and moments. 

An expression for the equivalent force F^ is 


F E = T + H 


(Cl) 


Since T and H are orthogonal, the magnitude of F^ is 


F c = | F_ I • (T 2 + H 2 )‘ /2 


(C2) 


by 


The magnitude of the total moment supplied by T and H is given 


M = M 


= (T 2 R t 2 + h 2 r 2) 17 2 

I n 


(C3) 


On the other hand, the magnitude of the moment due to F^ acting at 
distance R^ is 


" = «E F E 


Equating the two moment expressions, 


r e f e = (t2r t + H ' r h } 


1/2 


or 


R E - (T 2 R t 2 + H 2 R H 2 ) 1/2 / F e 
U sing equation (C2), 

R e = (T 2 R 2 + H 2 R 2 ) 1/2 / (T 2 + H 2 ) 1/2 


( C4 ) 


(C5) 


(C6) 


(C7) 
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or in terms of dimensionless quantities 


5 e = (T 2 ? 2 + h 2 ? 2 ) 1/2 / (T 2 + H 2 ) 1/2 

Equation (C8) represents the effective radius value used by the 
Subsonic Propeller Noise Module. 
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Figure 3.- Blade surface coordinate system. 


11.1-39 




Figure Al.- Roots of the retarded time equation. 
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11.2 TRANSONIC PROPELLER NOISE MODULE 


Sharon L. Padula 
Langley Research Center 


INTRODUCTION 

The Transonic Propeller Noise Module computes the periodic 
acoustic pressure signature and spectrum of a propeller with transonic 
tip speed. The computation is based on a solution of the Ffowcs 
Williams-Hawkings (FW-H) equation without the quadrupole source term. 
(See refs. 1 to 3.) The blade surface pressure is assumed to be steady 
and must be specified on a grid of points covering the blade surface. 
The observer is always assumed to be moving with the aircraft; only in 
such a frame is the acoustic pressure periodic. The solution to the 
FW-H equation is given by two different formulations that assume a 
source distribution on the full blade surface. The formulation that is 
suitable for blade sections moving at transonic speeds is used only 
when necessary. 


c 

f 

9 

k 

* 

£ 

M 

Nb 

Ns 

Nt 

n 

P 


SYMBOLS 

ambient speed of sound, m/s (ft/s) 

function defining blade surface S 

= 0 (describes relationship between source and 
observer time) 

harmonic number 

local force per unit area of blade acting on fluid. 
Pa (lbf/ft 2 ) 

source Mach number 

number of blades 

maximum harmonic number of interest 
number of time points in time history 
blade-surface normal vector 
surface pressure, p g - p^, Pa (lbf/ft 2 ) 
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p' 

Pl 

Pp 

P T 

Pa 

p 0 

Pref 


SPL 

t 

V F 

v 

X 

x 

y 

r 

n 

0 


p 0 


acoustic pressure. Pa (lbf/ft 2 ) 

acoustic pressure produced by loading. Pa (lbf/ft 2 ) 

total acoustic pressure produced by propeller. 

Pa (lbf/ft 2 ) 

acoustic pressure produced by thickness. 

Pa (lbf/ft 2 ) 

absolute surface pressure. Pa (lbf/ft 2 ) 

ambient pressure. Pa (lbf/ft 2 ) 

reference pressure for definition of decibels, 

2 x 10- 5 Pa 

distances from source point at emission time to 
observer, r = ix - y (x*)i, m (ft) 

sound pressure level with respect to p rg ^ dB 

time at which noise signal is received by 
observer, s 

forward velocity of aircraft, m/s (ft/s) 
source velocity, m/s (ft/s) 

observer position in aircraft-fixed frame, m (ft) 

observer position in ground-fixed frame, m (ft) 

source position in ground-fixed frame, m (ft) 

locus of points associated with emission time t 

source position in blade-fixed frame, m (ft) 

angle between radiation vector and surface 
normal vector, rad 

pitch change angle at root, rad 

elliptic surface coordinates 

ambient density, kg/m 3 (slug/ft 3 ) 
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T 


time at which noise signal is emitted at source 
position, s 

t* solution to retarded time equation, g = 0, for a 

given observer time, s 

<(> solution to retarded time equation as solved by 

Newton's method, n( T - t) 

angle between ai rcraft-fixed Xj-axis and 
blade-fixed nj-axis, rad 

q angular velocity of blade, rad/s 

Subscripts: 

i component along ith coordinate axis 

(note: r^ implies summation convention) 

n component in direction of surface normal 

0 ambient condition 

r component in direction of radiation vector 

Notation: 

-v vector 

A unit vector 

. time derivative 

* complex conjugate; specific solution to retarded 

time equation 


INPUT 

The computation of transonic propeller noise requires definition 
of the blade shape characteristics, aerodynamic characteri sti cs, flow 
and operating parameters, and observer coordinates. The noise 
signature is calculated by using two different full-blade formulations, 
depending on the radiation Mach number of the blade section. The input 
is consistent with that for the Subsonic Propeller Noise Module, and 
the two modules will give identical results for subsonic cases. 


11.2-3 


x + 


User Parameters 


N 

N 

N 

M 

R 


b 

s 

t 

z 


0 

* 


0 

0 


Q 


5 1 
5 2 

( 5 1 »?2 ) 

03(^1 
$1 ,$ 2 


5 1 

5 2 
*l> 

Pti lt 5 2 ,<P=0) 


number of blades 
highest harmonic number desired 
number of time points in time history 
aircraft Mach number 

blade length measured from axis to tip, m (ft) 
pitch change angle at root, rad 
initial azimuth angle of first blade, rad 
angular velocity of blade, rad/s 

Blade Shape Table 
spanwise location, re R 
chordwise location, rad 
blade surface abscissa, re R 
blade surface ordinate, re R 
slopes of basis functions at nodes 

Aerodynamic Characteristics Table 
spanwise location, re R 
chordwise location, rad 
azimuth location, rad 
pressure loading, re Pq(R^) 2 

Observer Table 

initial observer position, m (ft) 
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Atmospheric Properties Table 

Pg ambient fluid density, kg/m 3 (slugs/ft 3 ) 

c local sound speed, m/s (ft/s) 


OUTPUT 

This module provides the complex pressure spectrum and pressure 
time history at a fixed observer location measured from the center of 
the hub. Complex pressures are given with the e“‘> ^ time harmonic 
convention, and all spectra are understood to be two-sided with 
Pp(x»-k) = Pp*(x,k). With this convention, the mean-squared pressure 
is 2ppPp* (where the square root of 2 PpPp* is Pp ms ) for each 
harmonic. 


Complex Pressure Spectrum 
k harmonic number 

fk frequency, Hz 

C k (x) complex acoustic pressure, re pg(Rn) 2 

SPL k (x) sound pressure level associated with C k (x) 


Pressure Time History 
t t i me , s 

Pp(x,t) acoustic pressure, re p rms 


METHOD 


Acoustic Formulation 


The equation used to calculate the acoustic pressure time 
signature of a high-speed propeller (see ref. 4), derived in the 
Subsonic Propeller Noise Module, is 


4irp'(x,t) 


1 9 f 

p 0 cv n + V 

dS + / 
T* f=0 

1 

L_ 

<=>* 

1 

c at f J Q 

r | 1 - M r| 

_ r2 1 - M r|J 


dS (1) 

r* 
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Equation (1) contains a potentially singular term, |l - M r |, in 
the denominator. As long as each point on the blade movis subsonically 
in the direction of the observer, the solution behaves nicely. In this 
case, equation (1) can be written in the following computationally 
efficient form: 


4-rrp.' (X,t) = 7 - / 


f=0 


£ i r i 


r(l - M r ) : 


dS + / 

r* f=0 


V - £ i M i 

L^ 2 u-M r ) 2 


dS 


and 


.1 f 

c f=0 


i ( rM . r . + cM - cM 2 ) 
r i i r 


r 2 (l 


M ) 3 

r ' 


dS 


(2a) 


4irp|(x,t) = / 

1 f=0 


p 0 v n ( rM i ^i + cM r " 


r 2 (l - M p ) 3 


dS 


(2b) 


p‘(x,t) = p,‘ (x,t) + p ‘ ( x ,t ) 


(2c) 


where p^ and pj denote the loading and thickness noise, 
respectively. The symbols M^, M^, £. , and are defined as follows: 


M i = v i /c 

(3) 

• 

M i = V i / c 

(4) 

£ i = P"i 

(5) 

• A A 

£. = pn i + pn i 

(6) 


where 

n i = (n X n) i (7) 


When the radiation Mach number approaches or exceeds unity for a 
given point on the blade at a given time, then the following equation 
is used in place of equations (2): 
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Pp(x,t) = p^x.t) + p|(x,t) (Be) 

where e is the angle between the blade normal n and the radiation 
vector r. The symbol r is used to denote the locus of points on the 
blade that are at a distance r from the observer at emission time x. 
In other words, r is the curve of the intersection of the surfaces 
f = 0 and g = x - t + (r/c) = 0. 

The Transonic Propeller Noise Module uses equations (2) and (8) 
to calculate p'(x,t). Equations (2) are used whenever possible because 
the computer implementation is simple and efficient. A thorough 
discussion of these and other acoustic formulas for calculating 
propeller noise is contained in reference 3. 


Frames of Reference 

In this section the coordinate frames used in the acoustic 
calculations will be discussed. Basically, there are two reference 
frames: the ground-fixed x-frame and the blade-fixed n-frame. The 

x-frame, which remains fixed to the undisturbed medium, is set up as 
fol 1 ows: 

The origin of this frame is at the propeller center at time 
t = 0. 

The propeller disk is in the XjX 2 -plane and the x 2 -axis 

coincides with the propeller axis, with the positive direction 
in the flight direction. 

The Xj-axis is assumed to be orthogonal to x 3 , upward and 
parallel to the plane of symmetry of the aircraft. 

The x 2 -axis is defined in such a way that the x-frame is 
right-handed. 

The observer position is always specified in this frame as will be 
explained shortly. 
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The blade-fixed frame is set up as follows: 

The origin coincides with the origin of the x-frame at t = 0. 

The n 3 -axis coincides with the x 3 -axis. 

The n 2 -axis lies on the pitch change axis of the blade. 

The m^-plane is normal to the propeller shaft. 

The nj-axis is defined to make the n-frame right-handed. 

Note that the blade mean surface does not lie in the n^-plane but 
makes an angle equal to the geometric pitch angle at each radial 
position. Note also that the n-frame rotates with angular velocity n. 

In the acoustic formulas presented in the previous section, many 
operations involving dot products or cross products of two vectors are 
required. When the products involve blade normals or surface 
pressures, it is convenient to transform the components of these 
vectors to the n-frame and then perform the operation. However, note 
that the source position y always refers to a source point in the 
x-frame with a known positioner! in the blade-fixed frame. The 
quantities r = x - y and r = r/r are calculated in the x-frame, and 
then the vector components are calculated in the n-frame. Note that 
for any variable that has a velocity term, such as M r = r-j 
and v n = v • n, the velocity is with respect to the x-frame. In 
other words. 


The components of this vector must then be calculated in the n-frame to 
get M r and v n . 

Another frame of reference that is used conceptually is the 
aircraft-fixed , nonrotating t-frame. This frame coincides with the 
ground-fixed x-frame at t = 0. The acoustic pressure p‘()T,t) for an 
observer in the ai rcraft-fi xed frame with the position vector it is 
found as follows: Since both t and t are known, the relative position 

of the ground-fixed frame and the ai rcraft-f i xed frame is known. From 
this, calculate the observer position x in the ground-fixed frame. 
Using this x and t, calculate p‘(x,t), which is the same as the 
acoustic pressure at position )T in the ai rcraft-fi xed frame and at 
time t. Thus, the )t- frame itself is never used in acoustic calcula- 
tions, even though the pressure signal that results is p'()T,t). 
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The transformation from one frame to the other at time t 
written in matrix notation as 



can be 


(10a) 



where (ni,o 2 ,n 3 ) is a point in the^rotating n-frame, (X x ,X 2 »X 3 ) is the 
correspondi ng point in the moving X-frame, (x 1 ,x 2 ,x 3 ), is the 
corresponding point in the ground-fixed frame, and ijj (where ij>=rzt) is 
the angle between the Xj- and m-axes. 


A final coordinate transformation is necessary to account for 
changes in the root pitch setting. This rotation through a small 
angle e n , is applied to the blade surface coordinates that are 
calculated by the Blade Shape Module. Thus, 


n ( % i > ? 2 ) 



( 11 ) 


Whenever the n-frame or pitch-axis coordinate system is mentioned 
herein, it is assumed that the previously given corrective 
transformation has occurred. Often the Blade Shape Module and 
Transonic Propeller Noise Module will be executed without change to the 
collective pitch setting. Thus, 0g = 0 and the last transformation is 
unnecessary. 


Computational Strategy 

Three elements that influence the computational strategy for noise 
prediction are numerical accuracy, computational efficiency, and 
conceptual simplicity. The trade-offs are never straightforward, but 
the simplest method that gives numerically stable answers in a 
reasonable amount of computer time is desired. A good balance is 
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particularly crucial to the Transonic Propeller Noise Module. 

Moreover, it is important that the transonic and subsonic propeller 
noise methods give identical results for propellers with subsonic Mach 
numbers. These objectives guide the computational strategy. 


The first step in computing an acoustic pressure time history is 
to discretize one full revolution of the blade into even time 
increments t-j and to represent the blade by a set of source points 
n(j). Even though the tip of the blade is moving supersonically, some 
portion of the blade is subsonic. It will be helpful to define a 
subset of source points n$ that move subsonical ly . That is. 


TV - 




n 2 ,n 3 


(,Qn?) 




< 



( 12 ) 


For each time point, the noise produced by the part of the blade 
represented by r| s can be calculated by using equations (2). The 
numerical evaluation of equations (2) is discussed in detail in the 
Subsonic Propeller Noise Module. The noise produced by the supersonic 
tip of the blade is calculated by either equations (2) or equations 
(8), depending on the value of M r for each source point. If M r 
exceeds some transitional value (usually 0.98) for any source point, 
then equations (8) are used for the entire supersonic tip at that 
time. In practice, equations (8) are used for less than one-fourth of 
the time points even in the worst case. 

Equations (8) are evaluated by dividing the blade tip into panels 
and by defining a set of possible emission times { t j } . Each 
integral /f[] dr dx is approximated by applying a one-dimensional 
trapezoidal rulfe for a fixed xj and then applying another 
one-dimensional trapezoidal rule as xj varies. The partial 
derivative with respect to source time t is approximated by using a 
central difference formula after the integral over the supersonic tip 
has been calculated at times (t^ + At) and (t-j - At). 

The execution time required to compute equations (8) is greatly 
influenced by the set of emission times xj that are chosen. If the 
span of all possible emission times is divided evenly, a very large 
number of times will be required. A better plan is to determine one, 
two, or three time intervals during which the r curve intersects the 
current panel. Each of these intervals can then be divided evenly. 

The time intervals themselves can be determined by solving the retarded 
time equation at the four corners of the panel. This gives at least 
4, and at most 12, times xj at which the r curve is touching a 
single point on the panel. The time interval [xT, x*+l] will contri- 
bute to the integral over the current panel if the r curve is inside 
the panel during that time. It is sufficient to test the location of 
the curve at the median time (xj + xj+l)/2. If the distance from the 
observer to the r curve is in the same range as the distance from the 
observer to each of the four corners of the panel, then the time 
interval will contribute. 
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There are several sources of error in this prediction method that 
the user must understand. The most important involves the geometric 
model of the blade. The blade position vectors ri and normal vectors 
n describe the blade shape very smoothly and accurately. Each panel 
of the supersonic tip is actually a curved surface. However, since the 
trapezoidal rule is used to approximate the integrals and since linear 
interpolation is used to calculate the integrands at points along 
the edges of the panel, computationally the curved panel is replaced by 
a flat plate. It is very important that the user supply a grid that 
has a large number of nodes in those regions where the blade is highly 
curved. One good way to accomplish this is to use the same chordwise 
grid that was used for the Blade Shape Module but add extra points near 
the leading edge. It may be necessary to add extra spanwise grid 
points as well. It is wise to have at least five spanwise grid points 
in the subsonic region and five in the supersonic tip. Finally, the 
user must take care to specify enough time points N^ to represent 
the pressure time history accurately. Unfortunately, the computation 
time will increase dramatically as the number of time points and grid 
points increases. 


The Retarded Time Equation 

As discussed in the Subsonic Propeller Noise Module, The retarded 
time equation can be written in the form 

f(<t>) = A(}> 2 + B<(> + C + cos(<j) + 0) = 0 (13) 

where only the roots (<j>/o) < 0 are of interest. Figure 1 is a typical 
graph of the retarded time equation that often results when the source 
points are moving at close to sonic speeds. It is not clear from the 
figure whether this graph has one or more negative roots. In order to 
solve numerically for the roots using Newton's method, the number and 
approximate location of the roots must be determined. 

Notice that equation (13) could be easily solved if the cosine 
term were absent. The first step in finding approximate roots is to 
set cos(<(> + D) = -1 and solve for the negative root . The location 
of the actual roots must be in the interval [if>g,0] because the cosine 
term cannot have a value smaller than -1. Next, approximate the cosine 
curve by a parabola that has the same minimum point and the same 
intercepts, and then solve again. The resulting approximate retarded 
time equation may have one, two, or no real roots in the interval 
U ,0]. If there is only one root to the approximate equation, then it 
is a good initial guess to the one negative root of equation (13). If 
there are no real roots to the approximate equation, then ^ will be a 
good initial guess to the one root of equation (13). If, however, 
there are two roots to the approximate equation, then the situation 
must be further clarified. 
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Figure 2 is an exploded view of the graph in figure 1. Note that 
there are three roots in this particular case. The best way to find 
them is to find first the local minimum and maximum points on the 
curve. That is, solve 

■^ = 2A<(> + B - si n( <j> + D) = 0 (14) 

Approximate solutions, <j>j and <f> 2 , to equation (14) are adequate. If 
f ( 4 , j_ ) and f(<j> 2 ) have opposite signs, then there is a root of f 
between them. Similarly, if f(<j>g) and f(^ 1 ) have opposite 
signs, then a root exists between u them. If it is determined that three 
roots exist, then ( 4>g + <h)/2, ( <f>i + 4 > 2 )/ 2 » and <j> 2 /2 will be very good 
initial guesses for use by Newton's method. 
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11.3 PROPELLER TRAILING-EDGE NOISE MODULE 


Donald S. Weir 
PRC Kentron, Inc. 


INTRODUCTION 


This module computes the trai 1 i ng-edge noise produced by a 
propeller based on the method developed by Schlinker and Amiet (ref. 

1). The propeller geometry and flow characteristics are provided by 
the Blade Shape Module, the Blade Section Boundary-Layer Module, and 
the Propeller Performance Module. The observer is assumed to be moving 
with the aircraft, and the propeller is modeled as a noncompact 
source. The module produces a 1/3-octave-band spectrum due to the 
broadband noise source at each observer position. 

SYMBOLS 


b 

C 

c 

CO 

D 

E* 

F 

f 

9r 


L 


M 





blade section semichord, m (ft) 
blade chord length, re R 
ambient speed of sound, m/s (ft/s) 
directivity function 
combination of Fresnel integrals 
spectrum function 
frequency, Hz 

acceleration due to gravity, m/s 2 
(ft/s 2 ) 

wave number, u>/U, m" 1 (ft -1 ) 
chordwise integral of surface loading 
turbulence correlation length, m (ft) 
Mach number 

Mach number of airfoil section relative 
to fluid 

reference Mach number 
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4 



m 

N 


5i 

<p 2 > 


R 


R 


a 

r 


4 

r 


S 


s 

T 


U 

X,Y,Z 

->■ 

X 


x.y.z 

->• 

X 


Mach number of airfoil section relative 
to fluid 

reference Mach number 

Mach number of airfoil section relative 
to observer 

molecular weight of air 

number of points 

number of blades 

number of spanwise stations 

mean-square acoustic pressure, Pa 2 
( 1 b 2 /ft 4 ) 

blade length, m (ft) 

gas constant for air, m 2 /K-s 2 (ft 2 /°R-s 2 ) 

magnitude of r 

position vector, m (ft) 

surface pressure spectrum, Pa 2 -s 
(1 b 2 -s/ft 4 ) 

Strouhal number, f , 6/U 

I / O 

span of airfoil segment, m (ft) 
standard sea level temperature, K (°R) 
time, s 

velocity, m/s (ft/s) 

■f 

coordinates of X 

nonrotating coordinate system fixed to 
propeller hub, m (ft) 

coordinates of x 

coordinate system fixed to blade trailing 
edge, m (ft) 
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B 

6 

6 * 




0 


ni , n2 > ^3 

->• 

n 


0 

9 


0 


0 


0 


T 


X 


u 

Ci 


a 

♦ 

'P 
S 2 

w 

0) 

Subscripts: 

b 

c 

e 

f 


compressibility factor, (1 - M 2 ) 

boundary-layer thickness, m (ft) 

boundary-1 ayer displacement thickness, m 
(ft) 

= u( x 2 + 6 2 z 2 ) ^ /ct 
coordinates of n 

rotating coordinate system centered on 
propeller hub, m (ft) 

= 7 + H (M - x/a) 

polar directivity angle, deg 

root pitch angle, rad 

blade twist angle, rad 

wave number, u/U 

= Mk x /B 2 

spanwise position, re R 
ambient density, kg/m 3 (slug/ft 3 ) 

= fx 2 + B 2 (y 2 + z 2 ) ] ^ 

azimuthal directivity angle, deg 

blade rotation angle, rad 

blade angular rotation frequency, rad/s 

angular frequency, rad/s 

Strouhal number, a)6*/U 

bl ade 

convection 
emission time 
due to M. R. Fink 
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h 

i 

o 

t 

1/3 

Notati on: 


high-frequency approximation 
low-frequency approximation 
observer 
trailing edge 
1/3-octave-band value 

normal ized 
unit vector 
modi fied 


INPUT 

The computation of propeller trai 1 i ng-edge noise requires 
knowledge of the blade shape characteristics, flow and operating 
characteristics, boundary-1 ayer thickness, and observer geometry. The 
blade shape characteristics can be computed by the Blade Shape Module, 
the boundary-layer thickness can be computed by the Blade Section 
Boundary-Layer Module, and the flow characteristics can be computed by 
the Propeller Performance Module. 



User Parameters 

turbulence convection Mach number, re M 

reference Mach number, 

number of blades 

highest harmonic number desired 

blade length, m (ft) 

aircraft altitude, m (ft) 

root pitch change angle, rad 

array of blade initial azimuth angles, rad 
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5i 

CUj) 


9 T (?l) 

T ilt ( 5i ) 

n 3t (c i ) 



Blade Shape Data 
spanwise station, re R 
blade chord length, re R 
blade twist angle, rad 
blade trai 1 i ng-edge abscissa, re R 
blade trai 1 ing-edge ordinate, re R 

Observer Geometry 

array of observer coordinates, m (ft) 





Propeller Performance Data 

5i 

spanwise station, re R 

M(Ci) 

local Mach number 

<*Ui) 

angle of attack, rad 


Boundary-Layer Data 


spanwise station, re R 

a 

angle of attack, rad 

M 

Mach number 

«*(?!. a,M) 

boundary-layer displacement thickness, 
m (ft) 


Atmospheric Data 

z 

altitude, re R^T^/nw 

3 P r 

C (z) 
00 

ambient speed of sound, re c r 

p„(z) 

ambient density, re p r 
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Independent Variable Arrays 
array of spanwise stations, re R 


f array of 1/3-octave-band center 

1 / 3 frequency, Hz 


OUTPUT 

This module produces broadband noise data due to the source. 


f 

1/ 3 


-> 


X 

<p 2 (x,f 



i/3 


Broadband Noise Data 
1/3-octave-band center frequencies, Hz 


observer location, m (ft) 
mean-square pressure spectrum, 

on « ^ 


re p c 

00 00 


METHOD 

The methodology of this module is based on the work of Schlinker 
and Amiet (ref. 1), who present a technique for predicting helicopter 
rotor trail ing-edge noise. Several refinements have been made for 
improved accuracy and compatibility with the propeller prediction 
problem to include 

1. A change of coordinate systems to maintain consistency with 

other propeller analysis modules 

2. A modification of the problem definition to extend its 

validity to the geometric near field 

The presentation of the methodology is divided into several 
parts. First, a description of the theoretical noise prediction for a 
two-dimensional airfoil is presented. Second, an empirical 
two-dimensional airfoil noise prediction method is presented for 
comparison. Then, a description of the various options for 
directivity and spectrum functions is presented. Finally, the 
two-dimensional trai 1 i ng-edge noise model is applied to predict 
propeller noise. 


Theoretical Noise Prediction for Two-Dimensional Airfoil 

For the coordinate system shown in figure 1, the mean-square 
acoustic pressure per unit frequency due to trai 1 i ng-edge noise is 
gi ven by 


C-^ 
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<p 2 (x, 0 ))> 


( 1 ) 


k MzbV 

id 



£y(w) s ( oj) 


This result was developed in previous theoretical investigations by 
Amiet (refs. 2, 3, and 4). It is based on modeling the turbulent flow 
past the trailing edge of an airfoil as a frozen surface pressure 
pattern which convects downstream at a Mach number M c past the 
trailing edge. The pressure jump must be forced to zero at the 
trailing edge to satisfy the Kutta condition; thus, an induced pressure 
field that propagates away from the airfoil is produced. The details 
of the derivation of equation (1) are presented in reference 1. 

Several terms in equation (1) require further explanation. The 
wave number k x is defined as 


(A) 

u 


( 2 ) 


and the length scale a is 

a = [x 2 + B 2 (y 2 + z 2 ) ] 

The chordwise integral of the surface loading 


M 2 


(3) 


where 


L is expressed as 


Ll 


i20 


{!-(!+ i )E* [2( A + U M + c 0 )]} 


'X + yM + 5 \ l / 2 


ux/o + L 


(1 + i )E* [2(ux/a + c 0 )] 


(4) 


0 = A + u (M - x/a) 

A = ui/U 


U = M k x /B“ 


1/2 


c 0 = y (X 2 + 3 2 z 2 ) /a 


and the overbar denotes normalized by the airfoil semi chord b. 
function E*[x] is the complex conjugate of the integral 


(5) 

( 6 ) 

(7) 

( 8 ) 
The 


E[x] = / e 

0 


it 


dt 


( 2 7rt ) 


ttz 


(9) 
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Alternately, E*[x] can be written as a combination of Fresnel 
integrals as 


E* [x 1 = Ci(x) - i Si (x) (10) 

where Ci(x) and Si ( x) are the Fresnel cosine and sine integrals, 
respectively. 


The surface pressure spectrum function S( w) 
the empirical relation 

S(u) = 2 x 10" 5 (| pj} 2 ) 2 F ( to) 


is approximated by 


(ID 


where the function 

F ( w) 

is defined as 

E ( to) = 

(1 + 

u + 0.217 u 2 + 0 . 00562 to 4 ) 1 

and the Strouhal number 

to is 


( 12 ) 


0 ) 



(13) 


This surface pressure spectrum function was obtained from the empirical 
data of Willmarth and Roos (ref. 5). The spanwise turbulent 
correlation length 
(ref. 6) as 


£y( w) is estimated from data obtained by Corcos 


£y(w) - 


2.11) /oj 
c 


(14) 


By collecting all terms in equation (1) that are a function of 
position, a directivity function can be defined as 


D( 9, <t>,w) 


zr p 2. 

= ( -) L 

a 2 1 


(15) 


where the retarded source radius r e is introduced for dimensional 
consistency and 0 and <j» are the polar and azimuthal directivity 
angles, respectively (fig. 1). Substituting equation (15) into 
equation (1) yields 

, + /k Mb\ 2 s i (u) 

<p 2 (x, w)> = 1 0 ( 0 , 4> » (jo) S( to) (16) 

V 2 it / 2r 2 


Note that the directivity function 0 is a function of frequency in 
addition to the two directivity angles. It will be shown later that 
when properly normalized, the frequency dependence can be removed. 
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Equation (16) is the primary result of the theoretical analysis of 
reference (1). Its directivity and spectrum characteristics are 
discussed in later sections. In addition, comparisons are made between 
the theoretical prediction and the empirical t railing-edge noise 
prediction method. 


Empirical Noise Prediction for Two-Dimensional Airfoil 

Schlinker and Amiet (ref. 1) also present a scaling law for the 
trailing-edge noise of a two-dimensional airfoil. From analysis of the 
theoretical approach and available experimental data, they conclude 
that the dependence of the mean-square pressure on flow parameters is 

< p 2> x m 5 — D (17) 

r. 2 


where IT is the normalized directivity function given by 
D( e,<t»,o>)/D( ir/2 , tt/ 2 9 a>) • 

The normalized spectrum function developed by Fink (ref. 7) for 
1/3-octave-band data 


- , 3 / 2 -<+ 
F f (S t ) = 0.613(10S t ) 4 [(10S t ) + 0.5] 


where the Strouhal number defined as 


S 


t 



(18) 


(19) 


accurately predicted the available spectrum data for clean airframe 
noise where trailing-edge noise was assumed to dominate. The resulting 
prediction equation is 

5.279 x 10- 7 M 5 AL_ D F, (20) 

f 

where the proportionality constant was obtained from experimental data. 


<p 2 > 


1/3 


p„ 2 c 4 


Directivity Functions for Two-Dimensional Airfoil 

Both prediction equations (16) and (20) require usage of a 
directivity function. The three options available for a directivity 
function are 
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1. High-frequency limit to equation (15) 

2. Low-frequency limit to equation (15) 

3. Modified directivity function proposed by Fink (ref. 7) 

First, some analysis of the geometry of the problem is required. 
Equation (1) was derived in terms of the present airfoil position x, 
that is, the position of the airfoil when the sound is received. As 
can be seen from figure 1, for an airfoil moving in the -x-di rection , 
the observer coordinate at the time of sound emission is 

r cos 0 = x - Ut (21) 

e e v ' 


where t e = r e / c OT is the propagation time from source to 
observer. Thus, the present source position is 


x = r g (M + cos 0) 


( 22 ) 


and the y and z observer coordinates remain unchanged by the 
airfoil motion. Expressing these in the spherical coordinates of 
figure 1 yields 


and 


y = r g sin 0 cos 4) 


z = r sin 0 sin d> 
e T 


(23) 


(24) 


Substituting equations (22) to (24) into equation (3) for a yields 


a = r (1 + M cos 0) 


(25) 


and equations (22) to (25) allow equation (15) to be written as 


D ( 0 , <j> , to) = 


sin 2 e sin 2 < 
(1 + M cos 0) 1 


(26) 


Equation (26) serves as a more useful form to evaluate the four 
directivity options. 


Full theoretical directivity . Substituting equation (4) for the 
chordwise integral of the surface loading into equation (26) yields 
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D ( 0 , <f> , oj) = 


sin 2 8 sin 2 <t> 1 

(1 + M cos 0) 4 0 2 


e 1 20 {1 - (1 + i)E* [2(x + iM + ? Q ) ]} 


x + um + ? o y / 2 

yx/ a + / 


(1 + i)E* [2( yx/a + 



(27) 


The normalization of equation (27) requires the evaluation of 
D( tt/2, tt/ 2,0)). Substituting 0 = tt/ 2 and <j> = tt/ 2 into equations (22) 
to (25) gives x = Mr e , y 0, z = r e , and a = r e . Then from 

equations (5) and (8), 0 = X and e 0 = y; therefore. 


D( tt/2, tt/2,0)) = 


1 _ 

X 2 



(1 


(1 + i )E* [2 X + 2 y(M + 1) ]} 


+ 


( x + y(M + 1) V /2 
\ U(M + 1) ) 


2 


(1 + i )E* [2 y(M + 1) ; 


(28) 


Dividing equation (27) by equation (28) yields, for the normalized 
di recti vity 


TT( e , 4> ,to) 


sin 2 e sin 2 $ x 2 

(1 + M cos e) 4 o 2 


e i20 {l - (1 + i)E*[2(T + uM + c o )]} + (- 


X + pM + V 2 
UX/<J + 


(1 + i )E* [2(yx/a + c Q )] 


e 1 2 X { 1 - (1 + i)E*[2x + 2p(M + 1)]} + [-- j( - ffV - fp -] 1/2 (l + i )E* [2u(M + 1)] 


(29) 


Since X, y, t 0 , and 0 are all proportional to u>, the only 
frequency dependency in equation (29) is in the arguments of the 
exponential and E* functions. 

Equation (29) is obviously awkward and time consuming to use for 
multiple calculations. Fortunately, simplifying assumptions can be 
made without significant loss in accuracy. 

High-frequency approximation .- For a case where (Jo >> 1, 
corresponding to a large chord length or high frequency, the arguments 
of E* in equation (29) become large. Since 

1 im E* [x ] = 1 ^ (30) 

x->-°° 
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equation (29) in the high-frequency limit becomes 


Bu(e,+) = 


sin 2 9 sin 2 <p 
(1 + M cos 0 ) 4 


T 2 


"A + uM + 5 o ' 
A + y(M + 1] 


u(M + 1) 
(yx/a) + c 


(31) 


Note that in the high-frequency limit, the directivity function is 
independent of frequency. 

Additional insight into the high-frequency approximation can be 
gained by examining Dh(e,u/2). Setting <f> = -tt/2 , substituting 
equations (5) to (8) and (22) to (25) into equation (31), and 
rearranging yield 


D h (e,7r/2) 


2 sin 2 (6/2) 

(1 + M cos 0) [1 + (M - M ) cose] 2 


The normalization factor Dj-, ( tt/2 , tt/2 ) , which is required for 
substituting equation (31) into equation (20), is 


D h ( tt/2, tt/ 2, to) 


7 2 


1 - (M - M c ) 


(32) 


(33) 


where 0^ ( 0 ,<t>) = %{ */2,tt/2,u)) 0^(6, <t>). 

Low-frequency approximation . - For a case where ub << 1, 
corresponding to a small chord length or small frequency, the value of 
E*[x] approaches zero. Therefore, in the low-frequency limit, 
equation (29) becomes 




tt.Um) ■ 6 - s -A n --^- 4 

(1 + M cos 0) 4 0 2 


(34) 


Again, the directivity function is independent of frequency. Following 
the same logic as for the high-frequency case. 


D (0,tt/2) = ^7-7 

(1 + M cos 0) 2 [1+ (M - M ) cos 0] 2 


(35) 


D^( tt/2, tt/2, u>) 


A 2 


(36) 
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The low-frequency approximation of the normalized directivity 
function is also shown in figure 2 for comparison. Direct comparison 
between equations (32) and (35) shows that 


D h (e,u/2) 

D £ (e,7r/2) 


1 + M cos e 
1 + cos 0 


(37) 


and that in the limit as the Mach number approaches zero, the two 
normalized directivity functions become identical at <j> = ir/2. In 
addition to the difference in directivity shape, the two directivity 
functions have different normalizations. The ratio of equations (33) 
and (36) is 


D, ( tt/2 , tt/2 , to) 1 - (M - M ) 

D 4 ( ^2,11/2, u ) = (38) 

Thus, if the turbulence convection Mach number is set equal to the 
free-stream Mach number, the peak overall sound pressure level in the 
high-frequency approximation is 1/M larger than the low-frequency 
approximation overall sound pressure level. 

Modified directivity by Fink . - In the prediction method for 
airframe noise developed by Fink (ref. 7), the normalized directivity 
function for wing trai 1 ing-edge noise is 

D f (e s ,<t>) = sin 2 ( 0 / 2 ) sin 2 <t> (39) 


This function was modified for the Doppler amplification of a point 
source with velocity M in the ANOPP Airframe Noise Module (ref. 8, 
sec. 8.8) to yield 

D,(8,*) = sin 2 (e/2) sin 2 ♦ ( 

(1 + M cos 0) 4 

which, since Df ( tt/ 2 , -rr/2 ) = 1/2, results in the normalized 
di recti vity 


D f (e,*) = 2 s1 " 2 < e/2) s,n2 ♦ (41) 

(1 + M cos 0) 4 

This function is similar in shape to the other approximate directivity 
functions shown in figure 2. It differs in shape from the 
high-frequency approximation by the ratio 


D h (6,-rr/2) 

D f (e,Tr/2) 


(1 + M cos 0) 3 
[1 + (M - M ) cos 0] ; 


(42) 
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Note that the Fink directivity function Df has no frequency 
dependence, whereas D, D^, and are all proportional to l/oi 2 
before normalization. Since the Fink directivity function was obtained 
independent of the Amiet theoretical analysis, it is not surprising 
that it is not compatible with equation (16). Therefore, equation (41) 
will only be applied to the empirical prediction given by equation 
( 20 ). 


Spectrum Functions for Two-Dimensional Airfoil 

The theoretical noise prediction equation (16) has several terms 
which are a function of frequency. If it is assumed that the 
high-frequency approximation is valid, substituting equations (2), 
(11), (14), (32), and (33) into equation (16) and collecting all 
frequency terms yield 


<p 2 (x,oi)> 



(1.330 x 10- 7 ) 



[1 - (M - M c )]0 h (9,*) 2 ^ i 

(43) 


Equation (43) describes a narrow-band spectrum. It is often convenient 
to examine the 1/3-octave-band spectrum for comparison to experimental 
data. The definition of the 1/3-octave-band spectrum function is 


F i / 3 (w) = / 


F(S) 


du 


df 


(44) 


where £ and u refer to the lower and upper limits of the band. 
Evaluating the integral, assuming that the mean value of the spectrum 
equals the center value, yields for the 1/3-octave-band spectrum 
functi on 


F l/3 (5) = 0. 23077F ( w) 


(45) 


where Si is understood to be a 1/3-octave-band center frequency 
value. Therefore, from the definition of F(S>), 


0.23077 

1 + u + 0.217 Si 2 + 0.005625 4 


(46) 


or, expressed in terms of the Strouhal number and 
converting the left-hand side from an Si spectrum to a frequency 
spectrum 
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0.23077 


( 47 ) 


1/31 1 + 0.785S t + 0.134S 2 + 0.00213S£ 

To convert from 5 to S^, it is assumed that 6* = 6/8, which is 
typical for a flat-plate boundary layer. From equation (47), a 
normalized spectrum F 1 / 3 (St) can be defined such that 


* F i/ 3 (V = 


= 1 


(48) 


where the summation is over all 1/3 octave bands. In addition, the 
original spectrum function was a two-sided spectrum for -°° < £j < °°. 
Converting to a one-sided spectrum, multiplying by two for both sides 
of the airfoil, and applying the normalization condition yield 


F 1/3<V 


0.0754 

1 + 0.785S t + 0.134S t 2 + 0.00213S£ 


(49) 


Therefore, for a 1/3-octave-band analysis, the F(io)/o> term should be 
replaced with • r 1 / 3 (St)> where 


F l/3 (S t ) = 12.243 F 1/3 (S t ) 


(50) 


and F l/3 (S t ) is given by equation (49). 

Schl inker and Amiet (ref. 1) made comparisons of the spectrum 
function F(5) with available experimental data for the surface 
pressure spectrum of an airfoil. They determined that equation (11), 
which is based on flat-plate surface pressure data, tended to 
underestimate the airfoil data by as much as 7.0 dB. They determined 
the difference between the flat plate and airfoil data. This 
difference function aF(w) is defined in table I. They recommended a 
modified surface pressure spectrum function F'(uj) as 


F 1 ( ai) - F( u) AF( m) ( 51 ) 

This can be converted to a normalized 1/3-octave-band spectrum, 
following the procedure of equations (44) to (47), with the result 

0.0215 aF(S.) 

p 1 (s ) = i 

1/3 1 + 0.785S + 0.134S 2 + 0.0022is£ 


Therefore, for the Schlinker and Amiet modified spectrum function, the 
term F(u))/ui should be replaced by F'^^St), where 
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f ; /3 <v - 42 - 923 f ; /3 <v 


(53) 


In summary, the three normalized 1/3-octave-band spectrum 
functions available for trai 1 i ng-edge noise prediction are 


1. F i / 3 ( St ) given by equation (47) based on flat-plate 
surface pressure spectrum data 


2. F'i/ 3 (S t ) given by equation (52) based on airfoil 

surface pressure data 


3. Ff(St) given by equation (18) obtained by Fink (ref. 7) 

All three spectrum functions may be used for the empirical 
trai 1 i ng-edge noise prediction method. Only the first two can be used 
for the theoretical method since the constant for normalization of the 
Fink spectrum function is unknown. A comparison of the three 
normalized spectrum functions is given in figure 3. 

The 1/3-octave-band mean-square pressure for the full theoretical 
directivity and the low-frequency approximation can be determined in 
the same manner. In general, the 1/3-octave-band form of the 
theoretical prediction equation (16) is 

MbV S£ (u>) 

— J D(0,*,o.) S( w) (54) 

tt / 2r 2 
e 

where conversion to a one-sided spectrum and both sides of the airfoil 
are taken into account. The frequency w is assumed to be a 
1/3-octave-band center frequency value. 


<p 2 (x,w)> 


1/3 


0.923i 


Comparison of Two-Dimensional Prediction Methods 


The theoretical prediction in the high-frequency approximation 
given by equation (43) can now be directly compared with the empirical 
prediction method. Assuming that the convection Mach number M c 
equals the free-stream Mach number M, the boundary-layer displacement 
thickness 6* equals 6/8, and the modified spectrum function of 
Schlinker and Amiet given by equations (52) and (53) is valid, equation 
(43) is written in 1/3-octave-band form as 



= (7.136 x 10- 7 ) -^M 5 D h (e,<|>) F' l/3 (S t ) (55) 

r e 
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Using the hi gh-f requency approximation for the normalized directivity 
function D, the empirical prediction equation (20) is written as 

<p 2 > . _ _ 

iLi = (5.279 x 10- 7 ) — M 5 D.(e,<t>) F-(Sj (56) 

„ z_ 4 „ 2 n rt 


Comparison of equations (55) and (56) shows that the theoretical 
method predicts an overall mean-square pressure about 1.3 dB higher 
than the empirical method and the spectrum shapes differ as shown in 
figure 3. It should be noted, however, that equation (54), not 
equation (55), represents the full theoretical model. 


Geometry for Propeller Application 


Application of the two-dimensional airfoil trai 1 i ng-edge noise 
prediction to a propeller is accomplished by modeling the propeller as 
the sum of individual airfoil segments. The noise contributed by each 
segment is summed at the observer to produce the noise signature at a 
fixed time. As time changes, the observed noise spectrum changes 
because the relative motion of the blade about the hub. The resulting 
time-dependent; noise spectra must be further analyzed to produce the 
complete noise signature. 


The geometry of the propeller application is shown in figure 4. 

The nonrotating coordinate system fixed to the propeller hub is denoted 
by ? = (X,Y,Z). The observer is fixed relative to the hub and its 
location denoted by r 0 = (X 0 ,Y 0 ,Z g) in the t coordinate 
system. The position of the trailing edge of a specific airfoil 
segment is denoted by rp. The components of the trail ing-edge 
location are provided by the Blade Shape Module in the rotating 
coordinate system ri» where rp = (n-| t U 1 )»f;i,n 3 t U 1 )) and z L 
is the spanwise blade coordinate. The blade position rp must be 
transformed twice to express it in the t coordinate system. The 
first transformation to the n‘ coordinate system accounts for the 
blade root pitch e 0 as 


->■ 

r 

b 


cos e„ 0 sin e 
o o 


0 1 0 


-sin e 0 
o 


cos e 

o 



(57) 


The second transformation removes the rotation of the rj' coordinate 
system, which results in 
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where ^ is the blade rotation angle defined in figure 4. 

The observer position relative to the airfoil segment is now 
defined in the ^-coordinate system at the time of sound emission as 



The airfoil trail i ng-edge noise predictions are expressed in terms of 
the x = (x y,z) coordinate system. Thus, r e must be transformed 
to this coordinate system. The inverse of the transformations (57) and 
|58) can be used to transform to the n-coordinate system. Data in the 
x-coordinate system is transformed to the n-coordinate system in the 
Blade Shape Module by 



where 0j is the blade twist angle. Thus, r e is expressed in 
the x coordinate system as 


COS 

9 I 

0 

sin 8y 

-i 

COS 

CD 

O 

0 

sin e 0 ' 

0 


1 

0 


0 


i 

0 

-si n 

0 

0 

cos 0 


-si n 

0 

0 

cos 0 

_ 

T 


tJ 



0 


0. 


cos i|i -sin ij; 

sin cos iJj 

0 0 
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Substituting equations (57) to (59) into equation (61) and evaluating 
the inverse of the matrices yield 



COS (0J + 

O 

CD 

0 

-sin 

( 9 T + V 


cos 


sin 

1 > 

0 


X 

O 

t 

r e Ui,lO = 

0 


1 


0 


-si n 


cos 

<l> 

0 


Y o 


sin ( 0j + 

e o ) 

0 

cos 

( e T + e o\ 


0 


0 


0 


1 

O 

M 


COS 0y 

0 

-sin 0y 


n lt (51) 

0 

1 

0 


5i 

sin 0y 

0 

COS 0J 


n 3t (?i) 


(62) 


Note that equation (62) is conveniently written in terms of input 
quantities. 

The position of the airfoil at the observer time is 

p(5i»*) = PgUif’I'g) + cjM.0,0) At (63) 

where the airfoil is assumed to have rectilinear motion after sound 
emission. Since At = ( ip - ij> e )/n and t = equation (63) can 

be written as 

rUi,'!') = P e (5i,^ e ) + cjM.0,0) ( ip - ^ e ) (64) 

The time of sound emission t e = ip e /n is found from the retarded 
time equation 

t g + -jr- = t (65) 

oo 

or 

nr 

’"e = * " -ET 1 (66) 

oo 

where r e is the magnitude of r e . Thus, for each observer time, 
equation (66) is solved for the emission time t e = ^ e /n. Then, 
the present airfoil position is determined from equation (64). The 
resulting value of r is used for the airfoil segment noise predictions. 
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In the propeller application, the directivity angles 0 and <j> 
are not related to x, y, and z in the same manner as equations (22) 
to (24). They are defined in terms of the coordinates 
Tg - ( Xg , yg» Zg ) as 



cos 

e - i 



1 

e 


cos 

T r sin 0 

e 

The present 

source 

coordinates are 

X = 

r (M + 
e 

COS 0) 

y = 

r sin 
e 

0 COS <j> 

z = 

r sin 
e 

0 sin $ 


(67) 

( 68 ) 

(69) 

(70) 

(71) 


Note that the motion of the airfoil relative to the observer does not 
equal the motion of the airfoil relative to the fluid. 


Computation of Acoustic Pressure 

The 1/3-octave-band mean-square acoustic pressure for each airfoil 
segment can be computed by either equation (20) or equation (54), 
depending on whether the user desires an empirical or theoretical 
prediction. The desired directivity and spectrum functions must also 
be selected. Then, the total mean-square pressure at a given time for 
a single blade is 


<P 2 (X, 


,S 


+ ;t)> . 

1 1/3 


N^i 

<p2 (X o’ S t ; Sl ,i » t ^ > l/ 3 


(72) 


where the normalization by p£c 4 is assumed. The corrected Strouhal 
number S^' has been introduce? because it is necessary to account 
for the Doppler frequency shift due to the relative motion of each 
blade segment to the observer. This frequency shift is, in general a 
function of Si and ip. 


11.3-20 


The general equation for the frequency shift of a source in a flow 

of Mach number Mf and relative source to observer motion is 

r 


to 

i 

to 


1 


+ 




A 

r 

e 


(73) 


For the case of the motion of the airfoil segment. 


-► 






C 3t p 

00 


and 


(74) 


M f = (-M,0,0) ( 75 ) 


where r e is the unit vector in the r e direction. Substituting 
equations (74) and (75) into (73) and noting that r P =(x ,y ,z ) 
yield e e ,J, e’ e 


0) 

I 

Cl) 


1 + 


WjL fe + ^e % . h 

c .\ r e 3 * r e H r e j 

- Mx /r 
e 


Therefore, the corrected Strouhal number is 


(76) 


Q_ 

( 3X 

e 

9 *e 


3Z e . \ 

c. 

1 + “ 

a , cos e + 


sin e cos <j> 

+ sm 8 sin +j 



1 

- M cos e 



(77) 


For multiple blades, each with an angle ^ relative to the 
initial blade, the contribution of each blade is summed by the relation 


N b 

< P 2 (X 0 .s t ; t» i/3 - J <P 2 (X 0 ,s; : t *£)> 


1/ 3 


(78) 


where the initial blade position ^ is zero. Since equation (72) is 
evaluated at values of ip given by 


* = r 1 (i ■ 1} 


(i = 1, 2, ..., N) 


(79) 
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Equation (78) can be evaluated for an even number of blades by simply 
time shifting the single blade solution. For an odd number of blades, 
interpolation is required. The resulting time-dependent spectrum 
<p 2 (X Q ,St;t) > l j 3 represents the noise signature of the 
propeller to trai 1 i ng-edge noise which is averaged by the relation 


(3t 0 .s t )> 


/ 3 


2tt 



<P* ()t , S.; t)> dt 
01 1/3 


( 80 ) 


to yield the mean-square acoustic pressure at the observer. 
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TABLE I.- MODIFIED SPECTRUM FUNCTION aF(5) 


> 


(A) 

AF 

0.034 

1.000 

0.040 

1.585 

0.050 

2.292 

0.060 

2.883 

0.070 

3.467 

0.080 

3.803 

0.090 

3.980 

0.100 

4.170 

0.200 

4.787 

0.300 

4.364 

0.400 

3.803 

0.500 

3.629 

0.600 

3.467 

0.700 

3.312 

0.800 

3.161 

0.900 

3.161 

1.000 

3.021 

2.000 

2.631 

3.000 

2.292 

4.000 

2.088 

5.000 

1.997 

6.000 

1.904 

7.000 

1.820 

8.000 

1.737 

9.000 

1.659 

10.000 

1.585 






0.4 
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Polar directi 







Figure 4.- Coordinate system for propeller 
trai 1 i ng-edge noise prediction. 
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12. TONE PROPAGATION EFFECTS 


12.1 BOUNDARY-LAYER PROPAGATION MODULE 


Gerry L. McAninch 
Langley Research Center 


INTRODUCTION 

The presence of a boundary layer, hydrodynamic or thermal, over a 
surface causes refraction and reflection of any acoustic wave incident 
upon the boundary layer. Consequently, the sound field at the surface 
differs from that which would be expected if the boundary layer was not 
present. The Boundary-Layer Propagation Module calculates the complex 
acoustic pressure at a plane surface of uniform admittance due to a 
specified acoustic field incident upon the boundary layer. 


c 

oo 

c* 

F(g) 

U 

i 


K 


k x 

k* 

y 

L 

m 

N 

n 

P 


Pr 

P* 


SYMBOLS 

sound speed in undisturbed free stream, m/s (ft/s) 
local sound speed, re c^ 

Fourier transform operator 
fundamental frequency, Hz 

=/7l 

integer, number of points for which initial 
data are given 

x-component of wave propagation vector, re 1/6 

y-component of wave propagation vector, re 1/6 

distance over which initial data are given, re 6 

integer exponent of 2 

integer, number of harmonics of interest 

harmonic number 

local pressure, N/m 2 (lb/ft 2 ) 

reference pressure, N/m 2 (lb/ft 2 ) 

total acoustic pressure, re p c 2 

' OO 00 
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p* 

Pi 

T 

oo 

T* 

U* 

u* 

v* 

w* 

x,y,z 

AX 

x*,y*,z* 

x ,y ,z 
0 ’■ 7 0 ’ 0 

X 1 

x 2 

y s 

a 

B 

Y 
S 

6* 

e 

v 

V 


incident acoustic pressure, re p c 2 

1 00 CO 

incident acoustic pressure evaluated at y* = 1, re p c 2 

oo oo 

undisturbed free stream temperature, K (°R) 
local temperature, re T m 
local flow velocity, re c^ 

x-component of complex acoustic velocity, re c^ 
y-component of complex acoustic velocity, re c^ 
z-component of complex acoustic velocity, re c^ 
Cartesian coordinates, m (ft) 
distance between data points, m (ft) 
nondimensional Cartesian coordinates, re 6 
observer coordinates, m (ft) 

lowest value of x for which data are provided, m (ft) 
highest value of x for which data are provided, m(ft) 
location of singular point 

distance from singular point at which numerical 
integration is reinstated 

surface admittance, re 1/p c 

00 oo 

ratio of specific heats 

boundary-1 ayer thickness, m (ft) 

complex acoustic perturbation density, re p^ 

radius of neighborhood around singular point where 
Frobenius solution is used 

Fourier transform, with respect to x, of u* 

Fourier transform, with respect to x, of v* 

Fourier transform, with respect to x, of incident normal 

acoustic velocity at top of boundary layer 

Fourier transform, with respect to x, of the reflected 
normal acoustic velocity at top of boundary layer 
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n. 

i 


Fourier transform, with respect to x, of total normal 
acoustic velocity at top of boundary layer 

Fourier transform, with respect to x, of p* 

Fourier transform, with respect to x, of incident 
acoustic pressure at top of boundary layer 

Fourier transform, with respect to x, of reflected 
acoustic pressure at top of boundary layer 

Fourier transform, with respect to x, of total acoustic 
pressure at top of boundary layer 

Fourier transform, with respect to x, of p* 

local density, kg/m 3 (slug/ft 3 ) 

undisturbed free-stream density, kg/m 3 (slug/ft 3 ) 

nondimensional local density, re p 

oo 

circular frequency, re 6/c 


INPUT 

Input to the module must specify the incident complex acoustic 
pressure at the top of the boundary layer, the acoustic admittance of 
the surface at the bottom of the boundary layer, the boundary-1 ayer 
thickness, the local temperature within the boundary layer, and the 
properties of the undisturbed free stream. 

For each observer, the observer coordinates must be specified. 
Further, the complex acoustic pressure, associated with the incident 
wave, must be provided at K = 2 m equally spaced points along the 
line (x,y,z) = (x,6,z 0 ), where z 0 is the observer's z-coordi nate, 
and 6 is the boundary-layer thickness. The lowest value of x, x x , 
at which the pressure is provided, and ax, the distance between 
points, must also be specified, as must be K, the number of data 
points for which initial data are given. The observers should be 
located at x () = x x + (n - 1) ax. 

A table of the local temperature as a function of y for 
0 < y <_ 6 and a table of the surface admittance e as a function of 
f i s 7equi red. 

Parameters required are 6, U OT , T«,, and pc, the boundary- 
layer thickness, undisturbed free-stream velocity, temperature, and 
density, respectively. The frequency content of the incident wave is 
specified by supplying f x , the fundamental source frequency, and N, 
the number of harmonics of interest. 
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Initial-Data Table 


f 


VV z o 


P i (x,f) 


f 

e 


y 

T(y) 


6 

U 

00 

T 

00 

Poo 

fl 

N 

x i 

AX 

K 


frequency, Hz 

observer coordinates, m (ft) 

incident complex acoustic pressure, N/m 2 
(lb/ft 2 ) 


Surface Admittance Table 
frequency, Hz 

surface admittance, re 1/p c 

OO 00 

Boundary-Layer Mean-State Profile 
distance above surface, m (ft) 
local temperature, K (°R) 

Parameters 

boundary-layer thickness, m (ft) 

undisturbed free-stream velocity, m/s (ft/s) 

undisturbed free-stream temperature, K (°R) 

undisturbed free-stream density, kg/m 3 (slug/ft 3 ) 

fundamental frequency, Hz 

number of harmonics 

lowest value of x-coordinate, m (ft) 

distance between data points, m (ft) 

number of data points required 
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OUTPUT 


Output of the module consists of the complex acoustic pressure and 
the mean-square acoustic pressure at the surface as a function of x 
and frequency f. 


Acoustic-Pressure Table 


f 

x 

P*(f.x) 

<p*(f,x)> 2 


frequency, Hz 

position on surface, m (ft) 
complex acoustic pressure, re p c 2 

00 oo 

mean-square acoustic pressure, re p 2 c 4 

00 00 


METHOD 


The complex acoustic pressure at the surface is to be determined 
from a knowledge of the incident pressure at the top of the boundary 
layer, p* and the complex surface admittance, 6(f). The problem is 
illustrated in figure 1. The equations governing the propagation of 
harmonic acoustic disturbances in a stratified moving medium are 
(refs. 1 and 2) 


-i w *u* + U* v* + c* 2 = 0 


3X* 

dy* 

3x* 


av* 

U* 

+ c* 2 -^1 = 0 



3x* 

ay* 



+ U* 96 * 

dp* . , 

+ — — V* + p^ 

, ( 3U* 

+ i* 

3X* 

dy* 

3x* 

ay 

iEl . ( 

:* 2 (-i to* 6* + U* 

3 6* 

, dp* 

3X* 

n 

3x* 

dy* 


(1) 

( 2 ) 

(3) 

(4) 


Equations (1) and (2) are the x- and y-components of the linearized 
momentum equation. The z-component is not required, as it is assumed 
that the problem is strictly two-dimensional; hence, w* and the 
derivatives of all quantities with respect to z* are identically 
zero. Further, the body-force terms and the viscous terms have been 
neglected. Equation (3) is the linearized mass conservation equation, 
while equation (4) is obtained by assuming that the propagation process 
is isentropic. As a result of neglecting the gravitational force, 
c* 2 = 1/p*, since c 2 = yP/ p, and since both y and P are constant. 
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In equations (1) through (4), U*, c*, and p* are the mean flow 

velocity, local sound speed, and density, respectively, and all are 
functions of y* alone. The acoustic perturbation quantities are u* 
and v*, the x and y perturbation velocity components, 6*, the 
density perturbation, and p*, the acoustic pressure. 

Further, all quantities have been nondimensional ized, the 
velocities with c m , the speed of sound in the undisturbed free 
stream, x* and y* with 6, the boundary-! ayer thickness, and the 
densities p* and 6* with p , while 

oo 

ui* = n2 irf i ( 6/c ) (5) 

n 1 °° 

where f x is the fundamental frequency, and n, the harmonic number, 
is such that 

1 < n <_ N (6) 

and N is the number of harmonics of interest. 

The problem statement is completed by specifying the incident 
pressure at y* = 1 

p*(x*, 1; oi*) = p*(x*, oi*) (7) 

where p*(x*, u>*) is the complex incident acoustic pressure at the 
top of the bouRdary layer, p, (x,f), in nondimensional form 
and the admittance of the surface at y* = 0 

v*(x*, 0; oi*) 

p* ( x* , 0; oi*) 

where a change of variables from f to oi* has been accomplished for 
the independent variable of the boundary aamittance. 

The acoustic density perturbation may be eliminated between 
equations (3) and (4) to yield 

-i to* p* + U*^- + (— + — ) = 0 (9) 

n 3x* 3x* 9y* 

Equations (1), (2), and (9) form a complete set for the unknown 
functions p*, u*, and v*, which is to be solved subject to the 
auxiliary conditions given in equations (7) and (8). 


= -B( 


oi*) 


( 8 ) 
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The first step in the solution procedure is to perform a Fourier 
transform, in the variable x*, on equations (1), (2), and (9), and on 
the auxiliary conditions of equations (7) and (8). The transform of g 
is defined as 

f* 2 _ik*x* 

F(g) = e x (g) dx* (10) 

X 2 - X jl X 1 

with x 2 > x 2 . 

The resulting system of equations is 


(oj* - U*k*)y + iv - k*c* 2 n = 0 (11) 

n x dy* x 

(a>* - U*k*)v + ic* 2 n' = 0 (12) 

(«* - U*k*)n - k*y + i V 1 = n (13) 

where 

n(y*;kj,wj) = F(u*) (14) 

v(y*;k*,aj*) = F(v*) (15) 

n(y*;k*,a>*) = F ( p* ) (16) 


and the prime notation indicates ordinary differentiation of the 
dependent variables with respect to y*.- 

Equations (11) through (13) are to be solved subject to the 
auxiliary conditions 


n ( 1 ; k J , oj* ) = n 1 (k*;oj*) 

(17) 

where 


n i ( k* ; to* ) = F[ Pl (x*;u*) ] 

(18) 

and 

v(0;k*,u)*) = - 3 ( w* ) n(0;k*,oj*) 

(19) 

The unknown function y can be eliminated between equations 
(13) to provide the system of equations 

(11) and 

ic* 2 n' + (w* - u*k*)v = 0 

(20) 
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i(i 


,* _ 


U*k*) v* 


ik* 

x 


dll* 

dy* 


v + 


[(< 


v* - 


U*k*) 2 - 


c* 2 k* 2 ]ir =0 (21) 


The initial data, p*(x*;io*), must be given at K = 2 n evenly 
spaced points over the interval 0 x : £ x* £ x 2 , with the first point at 
x* = Xj. Use of the transform, as defined in equation (10), implies 
that the problem is solved for the periodic extension of p*(x*;w*) 
from -oo < x < <*>. This function has the period L = x 2 - x 1# Note 
that L is dimensionless and gives the number of boundary-1 ayer 
thicknesses over which data are provided. After the transform, 
n 1 (k*; w*) will be known for K values of k* in the range 

£ k* £ — Z — LL (22) 

L X L 


These values of k* are equally spaced at increments 2 tt/L . The 
integration of equations (20) and (21) from y* = 0 to y* = 1 with 
the auxiliary conditions 


n(0;kj, 

v(0;k*. 


or 


) * 1 

(23) 

) = -bU*) 

(24) 


ensure that the initial condition given by equation (19) is satisfied. 
The initial -value problem given by equations (20), (21), (23) and (24) 
can be solved by numerical integration, as long as the equation 

U* = w*/k* (25) 


is not satisfied for any value of y* in the range 0 £ y* £ 1. For 
the purposes of the current discussion, this will be assumed to be the 
case. If, at y* = y s , equation (25) is satisfied, y$ is a 
singular point of the system of equations (20) and (21). The 
discussion of this case is provided in appendix A. 

The solution of equations (20) and (21), subject to the initial 
conditions given by equations (23) and (24), provides values of 
n(l;k*,oj*) = ilj. and v(l;k*,w*) = \£, the transform of total 

acoustic pressure p and total acoustic velocity v at the top of the 
boundary layer. It is necessary to decompose n, and \l. 
into components propagating in the positive and negative y-di rections. 
These may be denoted the reflected and the incident fields, respec- 
tively. Thus, 

n T = n. + n r (26) 

and 


12.1-8 


(27) 


V T = v i + v r 


The quantities n. and v. are related to each other through the 
equation 


v. = 


(co* - U*k*) 


y*=l 


while n and v are related through 
r r a 


(28) 


v = 

r 


c* 2 k*n 
y r 

U* - u*k*) 


y*=i 


In these equations, k* is given by 


k* = {[(oi* - U*(y*)k*) 2 /c* 2 ] - k* 2 } 


1/2 


(29) 


(30) 


Since equations (28) through (30) are required only at y* = 1 where 
c* = 1, they may be written as 


v i = 


v = 


- k *y ”1 


(ui* - U*k*) 
v n x 7 


k* n 
y r 

(to* - U*k*) 
v n x 7 


y*=l 


y*=l 


(31) 


(32) 


and: 


k* = {[to* - U*(l)k*] 2 - k* 2 } 


1/2 


y n ' x J x 
Equations (26), (27), (31), and (32) can be solved for n. to yield 


(33) 


n. 


n T k* - («* - U*k*) v T 


2k* 

y 


(34) 


y*=i 


Equation 34 gives the magnitude of the (to*, k*) component of the 
pressure field which is required to produce a Component at the surface 
of unit magnitude. The actual magnitude of this component at the top 
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of the boundary layer is known from the Fourier transform of the input 
pressure and is given in equation (17). The value of n(0;k*, w*) is 
therefore given by x n 


n(0;k*. 




y*=l 


(35) 


Once n(0; k*,w*) is found for all values of k*, the inverse 
transform determi nes n the time-harmonic components of x the pressure field 
at the surface. 


The sound speed at any point within the boundary layer, in 
nondimensional form, is found from 


c*(y*) = (36) 

where T* is the local temperature T(y) nondimensional ized by T . 

oo 

In this module, the velocity within the boundary layer is given by 
the Pohlhausen (ref. 3) profile 


U 

U*(y*) = ( 2y* - 2y* 3 + y* 4 ) ^ (37) 

oo 

The mean-square pressure is given by 

<p*> 2 = 2p*p* (38) 

where the bar indicates the complex conjugate. Finally, the sound 
pressure level (SPL) is given by 

SPL = 10 log(<p*> 2 ) + 20 log( Pco C OT 2 ) - 20 log(p ) (39) 

where 

P r = 2 x lO” 5 N/m 2 (40a) 

or 

p r = 4.1777 x 10- 7 Ib/ft 2 (40b) 
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APPENDIX A 


SINGULAR POINT CALCULATION 

The discussion in the section entitled "Method" is presented with 
the assumption that the quantity (to* - U*k*) is never zero in the 
region 0 _< y* _< 1. When this is no? the case, equations (20) and (21) 
cannot be integrated numerically across the point y* = y where 

K - u * k ;> y s * °- 


In order to transfer the solution across the singular point 
y* = Ys> ^ is necessary to expand the functions n and v in series 
about the point y* = y s . The determination of these series is, 
perhaps, easier if the unknown functions v and v' are eliminated 
from equations (20) and (21) to obtain the second-order ordinary 
differential equation 

c* 2 (to* - U*k*) n' ' + 2 [c* 2 k* + c* (to* - u*k*) In' 

X dy* x dy* n x 

+ (to* - U*k*) [ ( to* - U*k*) 2 - c* 2 k* 2 ] n = 0 (Al) 

m x n x x 

for the function n. The function v is then determined from the 
equati on 

-ic* 2 n‘ 

v - (A2) 

(to* - U*k*) 

As a first step in the solution procedure for equation (Al), the 
coefficients of this equation are expanded in a Taylor series about the 
point y* = y . Thus, 


c* 2 K - u * k ;> s .f, a oj 


2 [c* 2 — — k* + c* — ( a,* - U*k*) ] = z a, . x J 
dy* x dy* n x j=0 


w . 

(to* - U*k*) [(to* - U*k*) 2 - c* 2 k* 2 ] = i a 2 -T J 


where 


t = y* - y c 
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Equations (20) and (21) of the section entitled "Method" are 
integrated numerically from y* = 0 to y* = y - e, 0 < e « 1. The 
solutions n(y s - e) and v(y s - e) are used to provide initial 
values for equation (Al). Here, 


iU* - U*k*) 

n'(y-e)= — " x 


-* 2 


|y* = y s -e 


(A7 ) 


The series solution of equation (Al) about y* = y s is then used to 
determine u(y s + a), while v(y s + a) is found from 



+ a) = 


- iC * 2 !! 1 

U* - U*k*) 


y*=y s + a 


(A8) 


Equations (20) and (21) are then integrated numerically from 

y* = y s + a to y* = 1, subject to the initial conditions given by 

equation (A8) and to the known value of n at y* = y s + a. 

(See fig. 2.) 

The parameter a is such that 

0 _< a _< e (A9) 

and 

y s + ex <_ 1 (A10) 

If the inequality 

y $ + e <_ 1 (All) 


is satisfied, a is set equal to e. Otherwise a is given by 


a - 1 



(A12) 


Since e << 1, it is assumed that sufficient accuracy is maintained in 
the solution for n, even if the series of equations (A3) and (A4) are 
truncated after the term multiplying x 4 . Thus, for 0 < t _< e, it is 
assumed that n is given by solutions of the equation 


(a 


01 T + a 02 T 




+ 


(a^g a^^x + a^x + a^^x + ^ 14 ^ + 

^ a 21 T + a 22 t2 + a 23 t3 + a 24 r4 ^ n = ® 


(A13) 
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The coefficients a-jj of this equation are related to the functions 
U*(y*) and c*(y*) in appendix 8. These coefficients are also 
presented in appendix C for the special case where U*(y*) is given by 
the Pohlhausen solution (ref. 3) for flow over a flat plate 

U 

U* = ^ (2y* 2 - 2y* 3 + y* 4 ) (A14) 

00 

and c*(y*) = 1 . 

The solution of equation (A13) in the neighborhood of x = 0 can be 
determined by the method of Frobenius (refs. 4 and 5) and is of the 
form 


n( t) = A n^x) + B TI 2 (t) ( A15) 

where n^x) and n 2 (x) are independent solutions given by (refs. 5 
and 6) 


M x) = x 3 z b 
n=0 r 


n 2 ( x) - dn 1 (x) Ln(x) + 1 + z c x n 


n=l 


Here, the constant d is given by 

“ [ a 21 ( a ll + a 02 ^ + a 22 a 01 ^ 


d = 


3a 2 
^01 


(A16) 

(A17) 


(A18) 


and the constants b n and c n can be determined through use of 

the equations presented in appendixes D and E. It should be noted that 

this solution is invalid whenever 


01 




y*=y c 


= o 


Equation (A19) can be satisfied only when 


dU* 

dy* 



= 0 


(A19) 


(A20) 


For the Pohlhausen velocity profile, this implies that y s = 1. If 
equation (A19) is satisfied, the solution procedure is terminated. 
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The values of n(y s +a) and v(y s +a) can be determined from 
equation (A15) once the constants A and B are found. These 
constants are related to n(y s -e) and v(y s -e) by 


c* 2 ^ - iU* - U*k*) vn 2 

c* 2 (n 1 i^ - n^n 2 ) 

c* 2 tui^ - iU* - U*k*) vn x 
c* 2 ^!^ - n^) 


(A21) 


(A22) 


The quantities H( y + a ) and v (y<-+a) can determined from 
equations (A15) and (A2). 

The solution is now complete except for locating the singular 
points. These points are given by the solution of the equation 

(to* - U*k*) = 0 (A23) 

n * y $ 


For the Pohlhausen solution, this implies that 


IV 


a)" U 

= = (2y - 2y 3 + y 4 ) — 

y , * ' y s y s' c 

J s k* “ 


Thus, 


Letting 


c u* 

^ - 2 * s 3 + 2y s - — = 0 

S S U „ k X 


<t> = 


C 0)* 

o° n 

U k* 

°° X 


it is clear that 

0 < <b < 1 


since 


0 < y < i 


(A24) 


(A25) 


( A26 ) 


(A27) 

(A28) 
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For y and <> within this range, the solution of equation (A25) for 
y s maf be written as (ref. 7) 


1/2 1/2 

[2 - z + 2/(1 + z) j 


(A29) 



( A30) 
(A31) 

(A32) 


The solution given is valid as long as <f> and y are restricted to 
the regions given in equations (A27) and (A28). s 
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APPENDIX B 


CONSTANTS a-jj FOR GENERAL CASE J 


a oi = - x ;( c * 2 u 1 ) 


1 02 ’ T 


(c* 2 U 2 + 4c*C 1 U 1 ) 


a 03 = “ k x t — 6 ~~ + (C 1 2 + C * C ^ ^ + C * C ^ 


★ . ★ ★ 


*04 ' 


£[c* 2 U 4 + 80 * 0^4 + 12U 2 (c 1 2 + c*c 2 ) + 8U x ( 3c : c 2 + 0*03)] 


a 10 " 2a 01 


a n = 2 (c* 2 U 2 + 0*044^* 


a 12 = ( 3U*c*c* + U 3 *c* z )k* 


[-(30^2 + 0*03)14 + 3(ci 2 + c*c 2 )U 2 + 50*0*11* + c* 2 U*] 


13 = T - 


a 14 = y [ ~ ( 3c* 2 + 4c*c* + 0 * 04 ) 14 * + (3c 1 *c 2 * + c*c 3 *)U*] + 
k* 

[8(c 1 * 2 + c*c 2 *)U 3 * + 7c*c 1 *U4* + c* 2 U b * ] 


a 21 = " k x 2a 01 


a 22 ~ ' k x 2a 02 


a 23 ■ -f[ a 03 * k x* U l* 3 l 


a 24 ' - k x 2 f a 04 4 < U i*V>k;] 


4he functions U* and c* are defined by the equations 

U* = ^V_ and c* = for i = 1, 2, 3, 4, 5. 

1 dy* 1 1 dy* 1 

2 A11 functions are to be evaluated at y* = y . 
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APPENDIX C 


CONSTANTS 


^or these 


a.. FOR POHLHAUSEN FLOW AND CONSTANT 1 SOUND SPEED 
J 


a oi - - k ;< 2 - 6 * s 2 + 4 ^ s 3 > 


a 02 = - k x<^s ' 6y s> 7 


‘03 


U 

■k* — (4y ) 
x c 5 


a^v/i = 

-k* — 

04 

x c 

00 

a 10 = 

" 2a 01 

a ll = 

- 4a 02 

a 12 = 

" 6a 03 

a 13 = 

" 8a 04 

a 14 = 

0 

a 21 = 

- k x 2a 01 

a 22 = 

" k x 2a 02 

a 23 = 

- k x 2a 03 + a 

a 24 = 

3a 01 a 02 ‘ 


equations, y $ is defined by equation (A29). 
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APPENDIX D 


FORMULAS FOR CALCULATION OF THE CONSTANTS b r 


-3(2a 02 + a u ) 


{4 ( 3aQ2 + a n^i + [ 3 ( 2aQ3 + + a 2l-l^ 


-{5(4 a 02 + a n )b 2 + [4(3a 03 +a^+a^b^[3(2a 04 + a^) + a 22 


_ -{(3a 14 +a 23 ) + [4(3a Q4 + a 13 ) +a 22 Jb^ [5 ( 4a Q3 + a 


For n >_ 5, b n is calculated as follows: 


-{(n + 2) [(n + l)a 02 + a n ]b n _ 1 a- [(n + l)(n a Q3 + a 12 ) + a 21 Jb p _ 2 ] 

(n + 3)na Q1 

{n[(n - l)a Q4 + a^] + a 2 2^ b n-3 {[( n - l)a 14 + a 23 ] b n _ 4 + a 24 b n-5^ 

(n + 3)na n (n + 3)na Q1 
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APPENDIX E 


FORMULAS FOR CALCULATION OF THE CONSTATNS c p 1 2 


Cj = 0 


Cj = Arbitrary 


■ { 3 ( a H + 2 a Q2 )c 3 + c 2 [2(a 12 + a Q3 ) + a ?1 ] + a g3 + d[S(a Q1 b 1 + a Q2 ) + a n ]} 


** { 4c 4 ( a 1 1 + 3a Q2^ + c 3 1 a l 2 + 2a Q3 ) + a 21-^ + c 2^ a 13 + a 04^ + a 22^ 


( d [ 7a oi b 2 + b l^ 7a 02 + a ll^ + 5a 03 + a 12 ] + a 2 


“ { 5c 5 ( a i i + % 2 ) + c 4 [ 4 ( a !2 + 3a Q3^ + a 2jJ * c 3t 3 ^ a 13 + 2a Q4^ * a 22^ + c 2^ 2a 14 + a 23^ 


d[9 a 0 i b 3 + (9a Q2 + a 11 )b 2 + (7a Q3 + a 12 )b 1 + 5a Q4 + a 13 l 


For n _> 7, c n is given by 


-(n - l)[ a n + (n - 2) 3g 2 ] c n _ 3 _ {(n - 2) [a 32 + (n - 3 ) a Q 3 ] + a 2l^ c n-2 


n(n - 3) a r 


n(n - 3)a r 


(Equation continued on next page.) 


3 The arbitrary constant c^ may be set equal to zero. 
2 The constant d is given in equation (A18). 
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{(n - 

3) [a 13 

+ (n - 


f a 22^ c n-3 



n(n 

- 3)a Q1 


{ [ ( n 

- a 14 + a 23 

] C n-4 + 

a 24 c n-5 1 



n(n - 

3) a Q1 


a f b n- 

. 3 ' 2 " ' 

3)a 01 + 

b n -4t< 2 " 

- 3 > a 02 + a lll + W (2n 





n(n - 3)a Q1 

d iV 


■ 7 >»04 

+ a 13 ] + 

b n-7 a 14 ^ 


n(n - 3)a Q1 
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12.2 TONE PROPAGATION MODULE 


Gerry L. McAninch 
Langley Research Center 


INTRODUCTION 

The transformation of pure tone data from the source frame of 
reference to the observer frame of reference and the determination of 
the sound field at each observer location is accomplished in a 
different manner than the same transformation of broadband data. The 
Tone Propagation Module performs this transformation for pure tone 
data, whereas the Propagation Module (PRO) must be used for the 
transformation of broadband data. 


The five steps required to transform pure tone data from the 
source reference frame to the observer reference frame, in the order in 
which they are performed, are 


1. Interpolation of the input source noise data, given as a 

function of emission time and observer directivity angle, to 
provide the mean-square acoustic pressure at the observer 
location as a function of observer time and source frequency 

2. Doppler shift of the source frequencies to give the 

mean-square acoustic pressure at the observer location as a 
function of observer time and observed frequency 


3. 


Determi nation of 
characteristic 


spherical spreading losses and 
impedance change effects 


4. Determination of atmospheric absorption losses 


5. Determination of ground reflection and attenuation losses 


The resulting output 
pressure as a function of 


table provides the mean-square acoustic 
frequency, reception time, and observer. 


SYMBOLS 

a incoherence constant 

C coherence coefficient 


c speed of sound, m/s (ft/s) 

f frequency, Hz 
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G 

H 

h 

k 

M 

M a 

n 


o 


<p 2 > 



<p 2 >* 


gr 


R 

R 

r 

Ar 

r 2 

SPL 

T 

t 

U(s) 

W(z) 

y 

a 

y 


ground effects factor 

acceleration due to gravity, m/s 2 (ft/s 2 ) 

altitude, m (ft) 

observer height, m (ft) 

wave number, 2-rrf/c, 1/m (1/ft) 

source Mach number at emission time 

molecular weight of air 

frequency harmonic number 

observer index 

mean-square acoustic pressure, N 2 /m 4 (lb 2 /ft 4 ) 

mean-square acoustic pressure accounting for 
atmospheric absorption 

mean-square acoustic pressure accounting for 
ground effects 

magnitude of spherical wave reflection coefficient 

universal gas constant, m 2 /(K-s 2 ) (ft 2 /(°R-s 2 )) 

distance, m (ft) 

path length difference, m (ft) 

source-to-image distance, m (ft) 

sound pressure level, dB 

temperature, K (°R) 

time, s 

unit step function 
complex error function 
coordi nate 

argument of spherical wave reflection coefficient 
elevation angle 
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n 

0 

0 

y 

y 

v 

e 

w 

p 

a 

♦ 

R 

Subscripts: 

a 

e 

o 

r 

s 

Superscri pt : 
★ 


dimensionless frequency; 2-npf Q /a 

incidence angle, deg 

polar directivity angle, deg 

absorption coefficient, nepers 

average absorption coefficient, nepers 

normalized ground admittance 

angle between source velocity vector and vector 
joining source and observer at emission time 

density, kg/m 3 (slug/ft 3 ) 

specific flow resistance of ground, kg/(m 3 -s) 
(slug/ (ft 3 -s ) ) 

azimuthal directivity angle, deg 
solid angle, sr 

ambient 
emi ssi on 
observed 

standard sea level 
source 

dimensionless quantity 


INPUT 

The input to this module consists of one or more source noise data 
tables computed for the same source coordinate system on the aircraft. 
These data tables consist of pure tone data. When two or more data 
tables are input they are merged prior to the transformation to the 
observer frame of reference. 

r $ source radius, m (ft) 

o specific flow resistance of ground, 

kg/(m 3 -s) (slug/(ft 3 -s)) 
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Af 


frequency bandwidth, Hz 


Source Noise Data Table 

f source frequency, Hz 

s 

0 polar directivity angle, deg 

e 

d, azimuthal directivity angle, deg 

t emission time, s 

e 

<p 2 (f g » 0 e »'t ) e * t e )>* mean-square acoustic pressure, re p 2 c^ 


t 

o 

f* e (t,o) 

t e (t,o) 

0 e (t,o) 

Y e (t,o) 

h(o) 

e w (t e ,°) 


y 

c*(y) 

p*c*(y) 

p*(y) 


Source Geometry Table 
reception time, s 
observer index 
distance, m (ft) 
emission time, s 
polar directivity angle, deg 
elevation angle, deg 
observer height, m (ft) 

Wind Axes Geometry Table 

angle between source velocity vector and vector 
joining source and observer 

Atmospheric Properties Table 
altitude, M a g p (H - H i)/ RT r 
speed of sound, re c p 
characteristic impedance, re p p c r 
density, re p p 
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Absorption Coefficient Table 
f frequency, Hz 

y altitude, M a g r (H - H 1 )/R'T r 

U average absorption coefficient, nepers/wavelength 

Engine Variable Table 
t source time, s 

M(t) aircraft Mach number 


OUTPUT 

This module produces a table of the mean-square acoustic pressure 
as a function of observed frequency, time, and observer. 


n 

t 

o 

c I<°> 

p*(o) 
f 0 (n,t,o) 
<p 2 (n,t ,o)>* 


Received Noise Data Table 
frequency harmonic number 
observation time, s 
observer index 

speed of sound at observer, re c r 
density at observer, re 
observed frequency, Hz 
mean-square acoustic pressure, re p 2 c 4 


METHOD 

Noise Data Interpolation 

Data tables from two or more sources in the same reference frame 
may be input simultaneously. If two or more data tables are input they 
are first merged, with components at equal frequencies summed, to 
provide a single table of mean-square acoustic pressure as a function 
of source frequency f s , emission time t e , polar directivity 
angle e, and azimuthal directivity angle <f>. That is, this input 
table provides <p 2 (f s , 0 e ,<j> e ,t e )>*. 

The transformation to the observer location is accomplished by 
using 0 e (t,o), <j> e (t,o), and t e (t,o), as provided by the 

Geometry Module, to convert <p 2 (f s ,0 e ,<j> e ,t e )>* to 
<P 2 ( f s,t,o)>*. 
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Doppler Shift 


Because of the relative velocity between the source and the 
observer, the observed frequencies are not the same as the source 
frequencies. The observer frequency is related to the source frequency 
through the equation 


y n > 

l - m cos e 


where M is the source Mach number, and e w is the angle the source 
velocity vector makes with the vector joining the source and observer, 
both quantities defined at the emission time, as shown in figure 1. 
This step yields a table of mean-square acoustic pressure, uncorrected 
for the various attenuation mechanisms, as a function of harmonic 
number and observation time for each observer <p 2 (n,t,o)>*. 


Spherical Spreading and Characteristic Impedance Change Effects 

The next step is to account for the attenuation due to spherical 
spreading and the change in mean-square pressure attributable to the 
difference in the characteristic impedance between the source and 
observer locations. These are determined by requiring conservation of 
acoustic energy. 

Figure 2 provides a schematic diagram of a conical ray tube of 
solid angle dfi. The acoustic energy passing through any cross section 
of the ray tube, per unit time, is the product of the acoustic 
intensity and the cross-sectional area of the tube. Conservation of 
energy considerations implies that this product must be independent of 
the cross section of the tube considered. This is expressed by the 
equation 

r 2 dn = 
s 


where r s and r 0 indicate source and observer locations, 
respecti vely . 

Equation (2) may be solved for the mean-square acoustic pressure 
at the observer location in terms of the mean-square acoustic pressure 
at the source location, the characteristic impedance ratio 
( pc) r s /£pc)r 0 » anc * the spherical spreading or area ratio term 

(r s /r 0 ) • 

In dimensionless form, this relationship is 

p*c*(y ) / r \ 2 

<p 2 (r o )> * = (r|) < " l(p s> > * (3) 
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where p*c*(y 0 ) and p*c*(y s ) are determined from the Atmospheric 
Properties Table for the observer height and source height, 
respectively. The dimensionless heights at the source and observer are 
given by 


y 


o 


= M g h/T R 
a 3 r r 


(4) 


and 


y $ = M a 9 r (r sin y + h)/T p R 


(5) 


The observer distance r, observer height h, and elevation angle y, 
are found in the Geometry Table. 


Atmospheric Absorption Losses 


The atmospheric absorption losses are now determined. The 
atmospheric absorption coefficient y is a function of frequency 
f 0 (n,t,o) and altitude y. The average absorption coefficient y 
is defined as 


u( f ,y s .y 0 ) 




_ y. u(f,y.) - y n u(f,yJ 

«y = (6) 


where_ y ? and y s are given by equations (4) and (5), respecti vely , 
and y is obtained from the Atmospheric Absorption Table. The 
equati on 


<p 2 >* 


abs 


= <p 2 >* 


-2(yf/c p )(r o -r s ) 


(7) 


provides the mean-square pressure at the observer. Equations (6) and 
(7) must be evaluated by using the observed frequency f 0 (n,t,o). 


Ground Effects 

The ground effect model used in ANOPP is based on the Chien-Soroka 
theory (ref. 1), coupled with the impedance function of Delany and 
Bazley (ref. 2). 


12.2-7 


The geometry required for the calculation of the ground effects is 
shown in figure 3. The source is located at an altitude H above a 
ground plane. The observer, located at a height h from the ground 
plane, receives sound which propagates along a direct path r and 
along a reflected path r 2 - The path length difference Ar = r^ - r 

is the most significant parameter used in the ground effects 
calculation. Use of the law of cosines shows that Ar is given by 

1/2 

Ar = (r 2 + 4h 2 + 4rh sin y) - r (8) 


The derivation of the Chien-Soroka theory is presented in 
reference 1. The resulting expression for the mean-square pressure 
with ground effects is 




[1 + R 2 + 2RC cos (a + k Ar) ] 


(9) 


where <P 2> *|abs 1S the free-field mean-square acoustic pressure, in 
this case corrected for spherical spreading, characteristic impedance 
change effects, and atmospheric absorption. Further, C is the 
coherence coefficient, and k is the wave number based on observed 
frequency. 


k = 2ir f o (n,t,o)/c(o) (10) 

where c(o) is the sound speed at the observer. The quantities R 
and a are the magnitude and argument, respectively, of the complex 
spherical wave reflection coefficient. The bracketed term in 
equation (9) 

G = 1 + R 2 + 2RC cos (a + k Ar) (11) 

is the ground effects factor. 

The coherence coefficient is the fraction of the initial acoustic 
energy which remains in phase throughout the propagation process. A 
reasonable approximation for the coherence coefficient is made by 
assuming a Gaussian distribution of the form 

C = exp r - ( ak Ar ) 2 ] (12) 

where a is the incoherence constant. Generally, a is given the 
value 0.01 which yields a value of 0.37 for C at k Ar = 32 tt. 
Substituting equation (12) into equation (11) yields 

G = 1 + R 2 + 2R exp [- ( ak Ar) 2 ] cos (a + k Ar) (13) 
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To account for a finite bandwidth, equation (13) must be averaged over 
the frequency width Af to provide the relation 

G . 1 + R 2 + 2R exp [- (ak An) 2 ] ( ak ar/2) (u) 

where Ak = 2n Af/c(o) . 

The complex spherical wave reflection coefficient, as determined 
in reference 1, is 


Re 1 a = r + (1 - r) F(t) 

where r is the plane wave reflection coefficient 

_ COS 0 - V 
COS 0 +V 

and F(t) is given by 


F( t) = 1 - Ar t W(i t) 


with 


and W 


1 / 2 

x = ( k r 2 / 2 i ) (cos 0 + r) 
is the complex error function 

co -t 2 

W(z) = y / * - t dt > °) 


(15) 

(16) 

(17) 

(18) 
(19) 


Further, cos 0 is given by 


cos 0 - r - s - V + 2h 
r 2 


( 20 ) 


funct 


For t 
ion al 


>10, an asymptotic approximation for the complex error 
ows F(t) to be expressed as 


F(t) - -2/tt U(-R e x) xe T + — 


2 x 2 


(2 x z ) 


2 ) 2 


( 21 ) 


where U is the unit step function 

r -\ 

1 (s > 0) 

U(s) = n/2 (s = 0)} 

0 (s < 0) 


( 22 ) 
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The remaining parameter to be determined is the normalized ground 
admittance v. The empirical expression 

- 0 * 75 - 0 * 73 -1 

v = [1 + (6.86n) + 1(4. 36n) ] (23) 

developed by Delany and Bazley (ref. 2) is used to determine v. In 
equation (23), n is the dimensionless frequency defined as 

n = 2 7 t pf q / o (24) 


Here p is the air density at the observer location and a is the 
specific flow resistance of the ground. A graph of the ground 
admittance is given in figure 4. 


With the relations for r and F(t), the magnitude and phase of 
the complex spherical wave reflection coefficient are computed through 
use of the equations 


R = | r + (1 - r) F(x) 
and 


(25) 


a = arg [r + (1 - r) F( t) ] 


( 26 ) 


For an acoustically hard surface ( n = 0), the theory is greatly 
simplified since n = 0, r = 1, R = 1, and a = 0. The expression 
for the ground effects factor (eq. (14)) reduces to 

a . 2 4. 2 exp[-(ak Ar)^] cos (k $ ^ (27) 


The mean-square pressure, with ground effects, is now calculated 
through use of the relation 


<p 2 >* 


gr 



(28) 


This final calculation provides a table of mean-square pressure at 
the observer locations as a function of harmonic number, observation 
time, and observer index, that is <p 2 (n,t,o)>*. 


User Options 

The Tone Propagation Module always performs the spherical 
spreading and characteristic impedance change calculations. The user 
has four options concerning the application of atmospheric attenuation 
and ground effects: 
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1. No atmospheric attenuation or ground effects 

2. Atmospheric attenuation only 

3. Ground effects only 

4. Both atmospheric attenuation and ground effects 

The standard output form for the noise data at the observer is the 
dimensionless mean-square acoustic pressure <p 2 >* as a function of 
observed frequency. 

Further, the user may request printed output of the sound pressure 
level (SPL), in decibels, defined as 

SPL = 1G log 10 <p 2 >* + 20 log 10 [p*(c*) 2 ] + 197 (29) 

ad 

Here p* and c* are determined from the Atmospheric Properties Table 
at the observer altitude y 0 . 


Applicability of Module 

In general the solution for the pressure field of a moving source 
contains two terms: one of which corresponds to the radiation field 

while the other is negligible at distances from the source which are 
large in comparison with both the source dimensions and the 
wavelength. The assumptions used to develop the theory used in this 
module imply that the observer is in the radiation field of the 
source. Further, it is assumed that the observer is much closer to the 
ground than to the source, that is h/r << 1. 
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Figure 4.- Empirical function for normalized ground admittance 




1. Report No. 2. Government Accession No. 

NASA TM-83199 , Part 3 

3. Recipient's Catalog No. 

4. Title and Subtitle 

Aircraft Noise Prediction Program Theoretical Manual - 
Propeller Aerodynamics and Noise 

5. Report Date 

June 1986 

6. Performing Organization Code 

505-31-33-11 

7. Author(s) 

William E. Zorumski and Donald S. Weir, Editors 

8. Performing Organization Report No. 

L- 1 5937 

10. Work Unit No. 

. .. . ..J 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, V A 23665-5225 

a * 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

14. Sponsoring Agency Code 


15. Supplementary Notes 


William E. Zorumski: Langley Research Center, Hampton, Virginia. 

Donald S. Weir: PRC Kentron, Inc., Hampton, Virginia. 

Part 1 (N82-19946) and Part 2 (N82-19947) of this document were published 
February 1982. 

16. Abstract 

Analytic propeller noise prediction involves a sequence of computations culminating 
in the application of acoustic equations. This paper describes the prediction 
sequence currently used by NASA in its Aircraft Noise Prediction program (ANOPP). 

The elements of the sequence are called program modules. The first group of modules 
analyzes the propeller geometry, the aerodynamics, including both potential and 
boundary-layer flow, the propeller performance, and the surface loading distribution. 
This group of modules is based entirely on aerodynamic strip theory. The next group 
of modules deals with the noise prediction, based on data from the first group. 
Deterministic predictions of periodic thickness and loading noise are made with time- 
domain methods. Broadband noise is predicted by a semi-empirical method. Near-field 
predictions of fuselage surface pressures include the effects of boundary-layer 
refraction and scattering. Far-field predictions include atmospheric and ground 
effects . 


17. Key Words (Suggested by Author(s)) 

18. Distribution Statement 

Aircraft noise 

Unclassified - Unlimited 

Noise prediction 
Propeller performance 
Propeller noise 

Subject Category 71 



For sale by the National Technical Information Service. Springfield. Virginia 22161 

NASA-Langley, 1986 













